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f ' 

What  began  as  an  evaluation  of  the  angular  integrals  arising  in 

) 

the  finite  element  solution  of  the  neutron  transport  equation,  has 

i 

grown  into  the  development  of  a computational  procedure  for  applying 
phase-space  finite  elements  to  the  same  equation.  The  entire  project 
is  an  extension  of  the  doctoral  dissertation  work  being  done  by  Capt 
John  Souders.  As  such,  I hope  that  this  research  answers  all  of  the 
questions  he  has  posed. 

I am  grateful  to  my  thesis  advisor  Lt  David  D.  Hardin,  PhD  for 
suggesting  this  topic  and  for  providing  advice,  support,  constant 
encouragement,  and  great  patience  throughout  this  endeavor. 

I am  most  appreciative  to  my  wife,  Paula,  and  my  daughters,  Chris-  i 

tine,  Amy,  and  Stacie  for  their  patience  and  understanding  in  enduring 
my  many  absences  while  this  thesis  was  being  prepared. 

Ronald  C.  Wheaton 
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Abstract 


Phase-Space  finite  elements  are  applied  to  the  static  rnonoenergetic 
neutron  transport  equation  in  two-dimensional  cylindrical  geometry  by 
computer  subroutines  written  by  the  author  to  collectively  assemble 
the  global  phase-space  matrix  for  solution.  The  technique  uses  a vari- 
ational formulation  based  on  the  second-order  self-adjoint  form  of  the 
transport  equation  within  which  the  dependent  variable  approximated  by 
the  finite  elements  is  the  even-parity  component  of  the  angular  flux. 
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APPLICATION  OF  PHASE-SPACE  FINITE  ELEMENTS  TO  THE 
NEUTRON  TRANSPORT  EQUATION  IN  CYLINDRICAL  GEOMETRY 

I Introduction 

The  label  finite  element  method  was  first  introduced  by  Clough  in 
his  treatment  of  plane  elasticity  problems  in  I960,  but  the  initial  ‘ 

efforts  to  use  the  method  had  appeared  even  earlier  in  the  applied 

mathematics  literature  with  the  work  of  Courant  (Ref  1:9).  Courant's  : 

approach  in  19^3  used  an  assemblage  of  triangular  domains  with  piecewise 

■ 

continuous  functions  defined  over  each  domain.  However,  his  novel 
approach  did  not  mature  until  Greenstadt  applied  it  and  developed  the 
mathematical  basis  for  the  finite  element  method  as  it  is  used  today 
(Ref  1:10). 

Unlike  finite  difference  methods,  which  lead  to  a pointwise  , * 

approximation  of  the  governing  differential  equations  over  a solution 
region  that  is  an  array  of  grid  points,  the  finite  element  method  uses 

V : 

a piecewise  approximation  to  the  governing  equations  for  a solution  > 

region  made  up  of  many  small,  interconnected  sub-regions,  or  elements.  j 

Over  each  element,  the  solution  of  the  governing  equations  is  assumed  \ 

to  be  both  a function  of  the  independent  variables  and  some  undetermined 
coefficients.  These  element  coefficients  are  then  determined  so  that 
the  assemblage  of  element  solutions  is  in  some  sense  an  optimal  approxi- 
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mation  to  the  true  solution.  In  addition,  the  coefficients  are  chosen  j 

so  that  any  required  continuity  of  the  solution  or  its  derivatives  is 

I I 

maintained  in  crossing  element  boundaries.  j 

' I 

By  the  early  1970’s  the  successful  application  of  the  finite 
element  method  to  problems  in  solid  mechanics,  heat  conduction,  and 
other  areas  led  to  the  application  of  the  method  to  neutronics  problems. 

Initially,  the  neutron  diffusion  equations  were  treated  most  extensively, 
followed  by  increased  app] ications  to  the  neutron  transport  equation  in 
one-dimensional  slab,  spherical,  ai.d  cylindrical  geometries  as  well  as 

two-dimensional  cartesian  geometry  (Refs  2;  3»  *0.  ‘ 

Solutions  of  the  transport  equation  based  on  finite  element  tech- 
niques  often  use  discrete  ordinal. o approximations  for  the  angular  ‘ 

variables  (Ref  J>).  Still  other  approaches  include  simultaneous  approxi- 
mations  of  both  the  angular  and  spatial  variables  by  finite  elements 
(Ref  2). 

The  purpose  of  this  paper  is  to  examine  several  aspects  of  the 
finite  element  method  as  applied  to  the  static  monoenergetic  transport 

equation  with  anisotropic  scattering  and  sources  in  two-dimensional  J 

cylindrical  geometry.  The  technique  uses  a variational  formulation 

based  on  the  second-order  self-adjoint  form  of  the  transport  equation 

within  which  the  dependent  variable  approximated  by  the  finite  elements 

is  the  even-parity  component  of  the  angular  flux.  This  component  is 

quite  attractive  for  finite  element  approximations  because  it  requires 

that  only  half  of  the  angular  domain  be  considered.  In  addition,  the 

even-parity  component  is  always  positive  and  easily  integrated  to  find 

the  scalar  flux  distribution  (Ref  1^9). 
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Before  discussing  the  computer  implementation  aspects  of  the  finite 
element  method  as  applied  to  the  monoenorgetic  neutron  transport  equation, 
the  theoretical  basis  for  the  second-order  even-parity  form  of  the  trans- 


port equation  will  be  developed.  Along  with  this  development,  a broad 
overview  of  the  finite  element  method  will  be  presented  and  the  vari- 
ational formulation  for  the  application  of  the  method  will  be  stated. 
Then  the  method  of  computer  application  will  be  discussed  by  describing 
the  computer  subroutines  written  by  the  author  to  assemble  the  system 
equations  for  solution.  Finally,  some  conclusions  will  be  presented 
regarding  needed  computer  code  improvements. 


3 


k 


II  Theory 


As  in  other  fields,  the  advances  in  computer  technology  have-  been 
the  driving  force  behind  the  scientist's  ability  to  solve  more  and  more 
complex  problems  as  well  as  the  means  by  which  he  can  increase  the 
accuracy  of  the  problems  he  has  already  solved.  But  advances  in  computer 
technology  are  not  the  only  means  of  increasing  accuracy.  Often  the 
application  of  another  numerical  method  can  bring  about  a significant 
increase  in  accuracy.  The  resulting  use  of  finite  elements  in  solving 
the  neutron  transport  equation  is  one  of  the  most  recent  steps  in  the 
quest  for  increased  accuracy. 

Present  discrete  ordinate  solution  methods  lead  to  anomalous 
scalar  flux  distortions  when  applied  to  transport  problems  having  strong 
absorbers  and  localized  neutron  sources  (Ref  5J?55)«  Thus,  as  more  and 
more  complex  geometries  are  encountered,  the  finite  difference  methods 
can  be  augmented  by  the  finite  element  method  whose  elements  can  accu- 
rately represent  complex  shapes  and  eliminate  the  ray  effects.  Conse- 
quently, computer  codes  based  entirely  on  the  finite  element  method 
have  been  developed  to  solve  the  transport  equation. 

In  this  chapter  the  neutron  streaming  term,  ib'V  (f  , of  the  trans- 
port equation  will  be  developed  for  the  problem  domain  being  considered 
in  this  thesis.  The  even-parity  form  of  the  transport  equation  will  then 
be  developed,  followed  by  a discussion  of  the  variational  formulation 
of  the  transport  equation  and  an  overview  of  the  finite  element  method. 

A brief  development  of  the  tensor  product  basis  will  conclude  the  chapter. 
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Problem  Domain  - Cylindrical  Coordinates 


The  monoenergetic  transport  equation  can  be  written  as 

^(r ,*)«  d> 

where  (f  is  the  streaming  term.  In  cartesian  geometry  this 

streaming  term  can  be  calculated  quite  simply.  It  has  the  form 

V\  - A*  '(cot,  + +^4l  (2) 

in  rectangular  coordinates  (Ref  6:59)  where  = SL ■ t and  X is  the 
angle  between  the  planes  formed  by  the  JL  and  z unit  vectors  and  the  z 
and  x unit  vectors.  However,  in  curvilinear  coordinates  the  streaming 
term  takes  on  a somewhat  more  complicated  form.  As  an  example,  the 
S--V<p  term  in  cylindrical  coordinates  is  (Ref  6:59) 

+^fL-^x(4|  - 4%)*/*  4*  (3) 

where  (f  is  the  polar  angle,  /*'*'*’  and  X is  the  angle  between  the 

planes  formed  by  the  JL  and  z unit  vectors  and  the  z and  r unit  vectors. 

Figure  1 shows  the  cylindrical  coordinate  system  corresponding  to  Eq  (3)* 

This  paper  deals  with  the  application  of  the  finite  element  method 

in  cylindrical  geometry  to  the  transport  equation.  The  problem  domain 

in  which  cylindrical  geometry  is  assumed  is  the  air  over  ground  burst 

problem  domain  (Ref  7)  of  weapon  physics.  Since,  in  this  domain  the 

air  density  varies  only  with  altitude,  azimuthal  symmetry  can  be  assumed 

and  the  conservation  form  of  Eq  (2)  becomes  (Ref  6:58) 

Vl- A1  C0&  X \ ^($y \-jk**iHTL) 

T *T  “ V ^90  * ^ Jfe  (4) 

Eq  (4)  is  the  streaming  term  which  will  be  calculated  in  this  thesis 

by  applying  the  finite  element  method. 
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Even-Parity  Transport  Equation  ( 

The  finite  element  method  will  he  applied  to  the  even-parity  form 
of  the  transport  equation  since  the  even-parity  component  of  the  angular 
flux  requires  that  only  half  of  the  angular  domain  be  considered,  making 
it  quite  attractive  for  finite  element  approximations.  In  addition,  the 
even-parity  component  is  always  positive  and  easily  integrated  to  find 
the  scalar  flux  distribution  (Kef  4:149).  The  following  derivation  is 
that  of  Kaplan  and  Davis  (Ref  8).  Consider  the  monoenergetic  transport 
equation, 

r,Jt)  + <£(v)^(r,  a)  = f ^(r,  £.£')%, -a'Mi'  + S(r„a)  (1) 

in  the  convex  region  D with  total  macroscopic  cross  section  and 

angular  scattering  cross  section  ( J:,  dt-Jp')  defined  such  that 

My)  (r, 

where  is  the  macroscopic  scattering  cross  section.  Assume  the 

region  is  surrounded  hy  a vacuum  so  that  on  the  surface  E of  D the 
vacuum  boundary  condition  is 

o , -s- -2«(t»)  < o (5) 

where  n is  the  outward  unit  normal  to  R and  rg  is  any  r £ R.  Eq  (l) 
is  also  valid  for  “A  , thus , 

Adding  Eqs  (l)  and  (6)  and  post -multiplying  by  1/2  leads  to 

+ i[S(r,a)  + S(  *:,-£)]  (7) 


i 


or 

with  the  terms  X(r  4$  * ± [^r.-fi)  - , ^4l)*  ♦ %,-i)], 

^(i,& ■*')«  i[os(x >-a-4/)  «-  0S(r , -&')  ] ' and  ^(x.-flr)  - [s(r  ,4*)  *S(crfl)] . 

Let  -*(*,*)  be  the  odd-parity  flux  and  n*,±)  be  the  even-parity  flux. 
By  substracting  Eq  (6)  from  Eq  (i)  and  multiplying  by  1/2,  we  get  the 
similar  result 

[<AM  + %,-•&)]  +^)  ^ [*M)  ‘ 

s / & fc(r>-  ■a0“  -&0] 

+ i[%&)-'s(rr^]  <9> 

or 

-& v ^,i)4  <s*(r)  X(r,4)*  / * &«£>*)  (10) 

where  the  additional  substitutions  STu^-a/)  * 

and  Mii»)  s-k[$fcr,i)~  %>--&)]  have  been  made.  Since  is  an  even 

function  of  &2.'  , that  is,  Cj(2C,-*--ar)  * , the  integral 

in  Eq  (8)  can  be  written  as 

+>faxw~-Vfe’'*y'*r/ 

= J i [^c,«04’ 

B / «j(r,4L-40  ^,*0 

and  Eq  (8)  becomes 

^c,*)  * / 4Sj(r,*)(11) 


; 


1 


I 


Similarly,  since  tSJ^  is  an  odd  function  of  QS\,'  such  that  4')  e 

- , Jt  -fl/)  1 the  integral  in  Eq  (10)  can  be  written  as 

= f *i,(x  *&■*')  "5  [%!*>-  <%£>-*:')]  dU' 

= / «i0c,*-*0  '*(*>*') 

and  Eq  (10)  becomes 

4-^  Hc,&)  + «*(r)  *(r,*)  * Jcr, *(r,  -Cr 12) 

By  putting  the  boundary  condition  in  terms  of  the  even-parity  and  odd- 
parity  fluxes,  we  can  write  Eq  (5)  as 

%.  > 4) : s H > *)  + ^(r. »-*)] + j [^(r*,-a)  - <%•>-■*)] 

- %*>*)  + ^r*,*) 

- 0 , -ft'fl  (I*)  < O (13a) 

where  Tg  is  any  r £ R and  n is  the  outward  unit  normal  to  R.  In  a 
similar  fashion,  we  can  write 

- V/(ri,*)-'X(ri,-ft )~o  , -&*(r,)>0  (13b) 

Eqs  (11),  (12),  (13a),  and  (13b)  are  the  first-order  form  of  the  trans- 
port equation  and  its  corresponding  boundary  condition  (Ref  8:l6?)« 

We  now  wish  to  determine  the  second-order  even-parity  transport 
equation  in  terms  of  ^ , the  even-parity  flux.  In  order  to  do  this,  we 
will  assume  Eq  (11)  to  be  an  integral  equation  for  and  introduce 

the  linear  operator  G (jr),  which  maps  functions  of  A into  other  functions 

6 

of  such  that  for  any  integrable  function, 
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Likewise,  introduce  the  G^(r)  operator  such  that 

GU2)  fc,4)]  * J<rUL(r,iv*/)4(r,4i/)A&/ 


(14a) 


(l4b) 


With  the  above  definitions  of  the  G^_  and  operators,  Eqs  (11)  and  (12) 
can  be  rewritten  as 


r Sj 


Gj1^  - Sj  - 4 vx 


(15a) 


4L  Su 


Gk%zSvl~ 


(15b) 


For  now  we  will  assume  that  the  G^  and  G^  operators  are  invertible  and 

-1  -i 

form  the  new  operators  K = G and  K = G . Operating  on  Eq  (15b) 

g g u u 

-1 

with  the  G^  operator,  we  form 

S;K*]  - 

or 

X * (16) 

Substituting  Eq  (16)  for  the  odd-parity  flux  into  Eq  (15a)  results  in 

S,f  . S,-  *-vfK„S*]  + 

or 

-*-v[Ka  (*'**)]  «■  <*  * * Sj  -*•*[!<*  Sa]  (17) 

which  is  the  second-order  even-parity  form  of  the  transport  equation. 
Its  corresponding  boundary  condition  is  the  dual  set 
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^,4)  4 KuJSu.  - = o , -£;-Tl  <0  ( l8a) 

kJs^  - &'V  f(x,£)]  = o , -n,-2L  >o  (l8b) 

I 

Earlier  we  assumed  that  the  G^  operator  was  invertible;  we  will 
now  determine  its  form  in  order  to  make  the  operator  physically 
significant.  The  operator  defined  by  Eq  (l4b)  can  be  proven  to  be 
3elf-adjoint  and  positive  definite  (Ref  8:174,175).  In  forming  its 
inverse,  we  will  take  advantage  of  the  fact  that  depends  only  on 

and  expand  .q.ju'J  in  Legendre  polynomials  as 

where  0^“-  and  aju  are  macroscopic  cross  sections  as  in  the  previous 

derivations.  Making  use  of  the  spherical  harmonics  addition  theorem  < 

(Ref  6:609),  we  can  write 

where  * denotes  the  complex  conjugate.  Substituting  this  expression 

into  Eq  (19)  gives,  ( 

00  4-£ 

and  Eq  (l4b)  expanded  in  spherical  harmonics  assumes  the  form 

* o 

- w*  ta  ' 

Collecting  terms,  we  have 

!«**•« 

11 


* 


or 


oo  +A 


G.®[W)]  !„,)  (20, 

where  S ^ are  the  Kronecker  delta  that  arise  from  the  fact 

that  the  Ylm  form  a complete  orthogonal  set  (Ref  6:608).  Eq  (20)  can, 


thus,  be  written  as 


G^f-ffc.4)]  Zj  Y^*)  - <*t  Y|J4.)f/w] 


or 


cy>  +1 


*1  £j<rt» -n'coR.Y^ 

The  inverse  of  G^(rO  is  the  integral  operator 

| qC  | 

K“[fM  i Gj  =E  lt  (<®«  - 0i"(r)j  ft.  r,.(u, 


(21) 


or 


a O 


idll 


*U?) 


0£(t)*j l«»mW  Cjtr) 

In  the  above,  the  cross  section  term  can  be  rewritten  as 


g;co- 


= I ♦ 


W-<£tn 


so  that 


12 


' . . •«% 

*; 


| 

n 

r*. 

t 

i 


* 

% 


or 


where 


Ku[<(t,«]=  ^ [ ka>  +/ ] 

. “ J£+t  t <Si^)  \ _ 

**  ( ^r)_  )f}(4*4 


(22) 


is  defined  as  the  macroscopic  scattering  cross  section.  It  will  have 
a finite  number  of  terms  whenever  the  scattering  cross  section  OjT  is 
expanded  in  a finite  number  of  terms.  A similar  expression  for  can 
also  be  derived  using  Eq  (l4a)  and  the  scattering  cross  section. 

Before  moving  to  the  section  on  the  finite  element  method,  we  will 
first  look  at  a variational  formulation  based  on  the  even-parity  trans- 
port equation  which  was  just  derived. 


Variational  Formulation 

As  has  beer  shown,  the  transport  equation  and  boundary  conditions 
can  be  derived  in  various  forms.  However,  a desirable  form  of  the 
problem  is  one  in  which  the  solution  minimizes  or  makes  stationary  a 
functional  subject  to  the  given  boundary  conditions. 

It  can  be  proven  that  the  solution  of  Eq  (17)  minimizes  the 
functional  (Ref  8:169) 


F(u.) 


Ktt(4'Vu)>  + <a*Gj  a> 
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-2<4t‘Vu,KuSa>-  2<“,S3>)  <*r jl  ol-a  j <li 

where  the  following  definition  of  the  inner  product  holds  on  the 
space  of  complex  functions  of^.: 

<-%0,W«)>  r /■ f to)  v>(a)  ju 


and  where  the  second  integral  is  a surface  integral  over  the  boundary 
fl  of  the  domain  D.  For  all  u satisfying  the  essential  boundary  con- 
ditions, the  variational  functional  results  in  the  following  equivalent 
weak- form  of  the  second-order  even-parity  flux  equation  as  derived  in 
Appendix  B 

<*1,G3  + >]  <Ll 

+ |jf  s J^f  <-&  •1^yl,KuSa> 

+ <Y1)S^>]  (24) 

The  finite  element  method  will  be  implemented  in  the  following 
chapters  using  this  equation  for  the  even-parity  flux. 


Finite  Element  Method 

The  numerical  technique  of  finite  differences  uses  difference 
operators  to  approximate  the  derivatives  in  a partial  differential 
equation.  In  contrast  to  this,  the  finite  element  method  does  not 
approximate  an  operator;  it  assumes,  for  an  assemblage  of  discrete 
elements,  a trial  solution  satisfying  the  boundary  conditions  of  the 
problem.  The  trial  solution  and  its  undetermined  coefficients  are 
required  to  satisfy  the  exact  equation  in  an  integral  sense;  this 
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determines  the  values  of  the  undetermined  coefficients. 

Additionally , the  accuracy  of  the  problem  solution  and  the  degree 
of  approximation  depend  upon  the  size  and  number  of  elements  used  and 
the  approximating  functions  selected  for  each  element.  The  following 
concepts  will  be  discussed  in  this  section:  the  element,  generalized 
coordinates  and  shape  functions,  natural  coordinates,  global  nodes,  and 
tent  functions. 

The  Element.  In  the  finite  element  method  the  problem  domain  is 
divided  into  a finite  number  of  subdomains,  or  cells,  which  are  inter- 
connected at  nodal  points  and  on  the  element  boundaries.  Of  course, 
this  division  can  be  quite  physical  in  nature,  with  each  cell  being 
thought  of  as  separate  from  another  like  building  blocks.  Or,  the 
division  can  be  mathematical,  with  the  problem  domain  (continuum)  being 
zoned  into  regions  by  imaginary  lines  or  planes.  No  matter  how  the 
division  is  done,  the  finite  element  method  solves  the  problem  collec- 
tively for  the  whole  domain  by  finding  a solution  for  each  of  its  parts, 
the  elements. 

Determining  the  shape  of  the  basic  element  to  be  used  in  the 
finite  element  method  depends  upon  such  things  as  the  problem  geometry 
and  the  number  of  independent  variables  needed  to  describe  the  problem. 
Thus,  one-,  two-,  and  three-dimensional  elements  with  straight  or 
curved  sides  are  possible.  Only  straight  sided  elements  such  as  those 
in  Figure  2 will  be  considered  in  this  thesis. 

As  stated  earlier,  the  element  .nodes  are  points  on  the  element 
boundary  where  adjacent  elements  are  considered  to  be  connected.  In 
addition,  they  are  the  points  where  the  field  variables  of  the  problem 
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are  used  to  define  the  approximating  functions,  also  known  as  interpola- 
tion functions,  for  the  element.  Therefore,  nodes  can  be  classified  as 
either  exterior  (boundary)  nodes  or  as  interior  nodes  (Figure  2c).  The 
element  interpolation  or  shape  functions  will  be  discussed  in  the  next 
section. 

Generalized  Coordinates  and  Shape  Functions.  In  the  finite 
element  method  those  func Lions  describing  the  behavior  of  the  field 
variables  within  an  element  are  called  approximating  functions,  inter- 
polation functions,  or  shape  functions  (Ref  1:131).  Polynomials,  which 
are  easily  differentiated  and  integrated  as  well  as  mathematically 
easier  to  handle  in  formulating  the  element  equations  and  in  computer 
calculations,  are  the  most  widely  used  shape  functions,  although  many 
other  functions  could  be  used.  Ultimately,  the  function  which  is  used 
should  obey  certain  inter-element  continuity  requirements  for  the  field 
variable  and  possibly  its  derivatives.  The  additional  requirement  that 
the  polynomial  expansion  remain  unchanged  during  a linear  transformation 
from  one  cartesian  coordinate  system  to  another  is  also  desirable.  In 
this  paper  only  two-dimensional  polynomials  will  be  used  to  generate 
the  shape  functions.  The  form  of  these  complete  two-dimensional  poly- 
nomials of  order  N can  be  written  as  (Ref  1:132) 

n.  . . 

j?*i 

where  the  number  of  terms  in  the  expansion  is  = (n+l)  (n+2)/2. 

As  an  example,  the  shape  functions  for  a rectangular  element  with 
nodes  at  each  of  its  corners  will  now  be  derived.  Consider  the  element 


in  Figure  3 over  which  an  expression  of  the  following  form  is  assumed: 

Cf* (x,  y)  = -*■  <*4  *y 


or 


where 


and 


r - M[«f 

M = hr  *1  *1  **4] 


(25) 


[p]  = [1 


v y xy 


] 


The  superscript  e indicates  that  the  equations  are  only  valid  for  the 
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e ‘ element.  The  coefficients  °^i  in  this  polynomial  series  representa- 
tion of  the  field  variable  are  called  generalized  coordinates  since  they 
have  no  physical  meaning  and  merely  fix  the  magnitude  of  the  solution  . 
An  evaluation  of  the  above  expressions  at  each  of  the  nodes  of  the  rectan- 
gle results  in  the  simultaneous  equations 


P\  * «■,'  + «2'V°'fh  + Vi 

tf\  - <X\  * ‘*2*2+'*3y2*  *4  *2y2 

<Pl  = <*f  + *2  *3*  “IV  *4  % 


or,  in  matrix  notation, 


where 


k‘]T  =[g*]  Mt 

{P}-  [ <P|  } 


M* 


i *,  y,  *»* 

I *z  ^z 

i 

I *4  >4  V4. 


The  generalized  coordinates  can  be  expressed  as  the  solution  of  Eq  (26), 
that  is, 


[«*]T«r6*jW 


Substitution  of  Eq  (27)  into  Eq  (25)  leads  to 


<P*  =[P][G*]'[(Pe]T-[N*][(Pe]T 


(28a) 


(28b) 


M =[P)[G*f 

Tlie  1L  are  the  shape  functions  associated  with  the  nodal  values  (or 
nodal  degrees  of  freedom).  As  a result,  the  undetermined  coefficients 
in  Eq  (28a)  now  have  physical  meaning  just  as  do  the  unknowns  in  a 
finite  difference  scheme.  Also,  the  shape  function  N.  referring  to 
node  i is  equal  to  one  at  node  i and  zero  at  all  of  the  other  nodes 
of  the  element. 

One  difficulty  encountered  in  calculating  shape  functions  in  this 
way  is  the  computational  effort  required  to  obtain  [ 6]  when  and  if 
it  exists.  Thus,  researchers  tried  to  obtain  the  shape  functions  by 
inspection  and,  as  we  shall  see,  they  succeeded  with  the  aid  of  a 
special  coordinate  system  called  natural  coordinates. 

Natural  Coordinates.  In  contrast  to  a global  coordinate  system 
which  is  defined  for  an  entire  body  or  structure,  one  can  define  a 
local  coordinate  system,  called  a natural  coordinate  system,  which 
applies  to  a specific  element  (Refs  1:138,  139;  9:88).  It  is  usually 
set  up  so  that  some  of  the  natural  coordinates  are  equal  to  one  at 
primary  external  nodes.  In  this  way,  the  natural  coordinates,  when 
used  to  derive  the  shape  functions,  not  only  simplify  the  formulation 
but  also  facilitate  the  evaluation  of  the  integrals  arising  in  the 
element  equations.  An  example  of  a natural  coordinate  system  for  a 
quadrilateral  element  (forming  a canonical  element  in  the  local  system) 
is  shown  in  Figure  4.  The  rectangular  element  of  Figure  4 is  also  a 
member  of  the  serendipity  family  of  elements  (Refs  1:170;  9:31~42) 
which  contain  only  exterior  nodes.  The  shape  functions  for  these 
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Fig.  h Natural  Coordinates  for  a General  Quadrilateral 

elements  were  found  by  inspection  and  have  the  form 

HLnY  T0+£;)(i+7Ph)  (29) 

for  a linear  element  (Ref  1:171).  The  element  maintains  the  necessary 
continuity  of  the  field  variable  (fi  along  its  boundaries. 

In  concluding  this  section,  it  should  be  noted  that  the  ease  with 
which  the  shape  functions  can  be  found  for  the  rectangular  elements 
makes  them  appealing  for  use  in  the  finite  element  method.  But,  their 
use  is  limited  because  they  cannot  represent  curved  boundaries  as  well 
as  triangles  or  elements  with  curved  sides. 

Global  Nodes  and  Tent  Functions.  Once  all  of  the  element  proper- 
ties have  been  found  for  a system  modeled  by  the  finite  element  method, 
the  overall  system  properties  can  be  determined  by  an  assembly  process. 
That  is,  the  element  matrix  equations,  Eq  (2tia),  describing  the  element 
properties  are  combined  to  form  the  matrix  equations  describing  the 
properties  of  the  entire  modeled  system.  The  assembly  process  uses  nodal 
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compatibility  as  its  basis;  the  value  of  the  field  variable  at  a node 
where  elements  are  connected  is  the  same  for  each  element  sharing  the 
node.  Figure  5 illustrates  this  point  for  a node  shared  by  five 
triangles.  In  order  to  implement  the  assembly  process  on  a digital 


computer,  a global  numbering  system  for  the  nodes  must  be  formed.  Figure 

6 shows  a mesh  of  four  linear  rectangles  along  with  its  global  node 

numbering  system  and  each  elements  local  node  numbering  system.  A 

global  node  number  I is  assigned  to  each  node  in  the  mesh.  In  each 

element  in  which  a given  node  appears  the  index  I is  used  to  identify 

it.  For  computer  use,  the  global  nodes  are  usually  stored  in  a matrix. 

Thus,  a matrix  IE  (e,i)  may  be  defined  such  that  IE  (e,i)  is  the  global 

th  fch. 

node  number  of  the  1 1 local  node  in  the  ec  element.  The  IE  matrix 


22 


Fig.  7 IE  Matrix  for  a Rectangular  Element  Mesh 
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and  the  global  node  values  of  the  field  variable  as  , then  the  form 
of  the  global  problem  solution  for  a mesh  of  II  global  nodes  is 

N 

(P{*,y)  * E <Pj  Tj(y,y)  (30) 

I*  I 

Tensor  Product  Basis 

In  this  section,  the  four-dimensional  space-angle  subspace,  within 
which  the  finite  element  method  will  be  applied,  is  formulated. 

In  the  work  of  Kaper,  Leaf  and  Lindetnan  (Ref  2:20),  the  finite 
element  method  is  used  to  solve  a six-group  transport  equation  for  the 
tensor  product  of  a finite  dimensional  subspace  whose  approximating 
functions  depend  on  the  spatial  variables  (x,y)  and  another  finite 
dimensional  subspace  whose  approximating  functions  depend  on  the  angular 
variables  Gm>$)  * In  this  thesis,  a similar  tensor  product  is  used, 
hut  the  spatial  and  angular  subspaces  are  generated  in  the  cylindrical 
geometry  previously  discussed.  Thus,  the  resulting  tensor  product 
subspace  is  four-dimensional  with  (r,z)  spatial  variables  and(yO.(X) 
angular  variables.  Mathematically,  the  formulation  can  be  written  as 
(Ref  10:135) 

N M 

X t Si(r,€)A:(/i,x) 

j*l  J 

where 

= functions  of  the  spatial  subspace 
Aj(M  = functions  of  the  angular  subspace 
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III  Computer  Application 

In  this  chapter  the  finite  element  method  will  he  implemented  in 
the  cylindrical  coordinate  system  of  chapter  II  using  the  weak  form  of 
the  even-parity  flux  equation 

JD  [<4-VT,Ku.(A.r*')>+  <>1,G}  f > ] Jr 

^ •/  [<Ji'V*t>KuSlt> 

+ > Si  ^ ) A*  (24) 

In  addition,  it  will  be  assumed  that  the  trial  functions  in  the  above 

formulation,  ^ , are  equal  to  tt  Vi  S;K»)  AAyU,x) 

i i*»  1 

and  that  the  test  functions,  >J  , are  equal  to  A^/x) SK(v, t)  where 
S is  the  spatial  shape  function  for  local  nodes  i and  K and  A is  the 
angular  shape  function  for  local  nodes  j and  L.  These  shape  functions 
are  calculated  by  use  of  Eq  (29).  By  choosing  the  test  functions  to  be 
the  same  as  the  trial  functions  used  to  represent  the  even-parity  flux, 
the  Galerkin  method  for  deriving  the  finite  element  equations  is  estab- 
lished (Ref  1:10tf). 

The  implementation  of  the  finite  element  method  in  Eq  (24)  is 
based  on  a modular  programming  approach  which  uses  FORTRAN  IV.  In 
this  approach  the  program  or  code  is  written  in  units  or  modules,  each 
one  performing  a basic  task  independent  of  the  other  modules.  In 
FORTRAN,  these  modules  are  subroutine  or  function  subprograms.  In 
addition,  the  finite  elements  used  in  the  formulation  of  the  element 
equations  are  the  4-dimensional  space-angle  phase-space  elements  of 
chapter  II  generated  over  the  tensor  product  subspace.  Since  it  is 
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difficult  to  visualize  a 4-dimensional  element,  the  subroutines  used  to 
code  the  even-parity  equation  deal  with  separate  rectangular  (r,z) 
spatial  and  angular  finite  elements  and  their  products.  These 

rectangular  space  and  angle  finite  elements  are  mapped  into  the  canonical 
rectangle  of  Figure  4. 

The  linear  transformation 


Vo 

1 

= 2 

y y 

r 

(Vr»)/(V4)  ‘ 

L'J 

L (vv.)  (v*,) 

l 

takes  a rectangle  having  an  arbitrary  set  of  coordinates  (r^z^), 
Cr-^Zp),  (r^,z-j),  (r/+  ,z^)  into  the  canonical  rectangle.  A similar 
transformation  is  used  for  the  (yUs Tt)  angular  element.  The  inverse 
mapping  is  given  by 


(v*i) 


131) 


for  the  spatial  element  and  a similar  form  is  used  for  the  angular 
element. 

Eq  (31)  and  a similar  equation  for  the  angular  variables  were 
implemented  on  the  computer  in  such  a manner  that  the  coordinate 
positions  in  the  global  spatial  and  angular  elements  could  be  deter- 
mined immediately  after  the  positions  in  the  corresponding  canonical 
elements  had  been  assigned.  This  logic  was  used  so  that  all  the 
element  computations  could  be  done  for  the  canonical  element  with 
changes  to  the  global  variables  made  as  needed.  As  implemented  on  the 
computer,  this  procedure  also  permits  calculation  of  the  spatial  and 


angular  shape  functions  at  local  nodes,  i,K,j,L  in  the  canonical 
elements  corresponding  to  the  global  space-angle  elements  being 
considered.  In  addition,  the  modular  nature  of  the  shape  function 
routine  allows  either  the  use  of  higher  order  elements  or  the  use 
of  entirely  different  elements  without  cfianging  any  of  the  other  sub- 
routines. 

As  stated  earlier,  the  coordinate  limits  of  the  canonical 
element  (-1  to  + 1)  make  Gaussian  numerical  integration  most  attrac- 
tive for  evaluating  the  space-angle  integrals  of  Eq  (24).  As  a 
result,  Gauss-Legend re  quadrature  was  used  in  the  subroutines.  As 
implemented,  the  numerical  integration  is  carried  out  on  an  element 
by  element  basis  over  the  four  space-angle  variables  simultaneously. 

The  integrands  calculated  in  the  procedure  are  the  inner  product  and 
boundary  terms  of  Eq  (24)  evaluated  for  all  possible  tensor  products 
of  the  local  shape  functions.  The  resulting  integrands  are  thus 
four-dimensional.  Since  CDC  FORTRAN  allows  the  use  of  no  more  than 
three  dimensions  in  an  array,  the  integrand  arrays  had  to  be  compacted 
into  two-dimensional  arrays  in  the  integration  subroutine. 

Provisions  have  also  been  included  in  the  subroutines  for  use  of 
any  combination  of  up  to  five  user  stored  Gaussian  quadrature  rules. 

In  addition  to  the  space-angle  integrations,  a second  integration 
over  the  entire  angular  domain  is  required  by  Eqs  (14a)  and  (22)  for 
the  G and  K operators. 

e u 

A closer  look  at  Eqs  (l4a)  and  (22)  reveals  that  for  every  space- 
angle  element  being  integrated  the  operators  each  introduce  an  addition- 
al integral  which  involves  not  only  the  current  neutron  direction 
but  also  all  other  possible  jj/  directions.  For  this  reason  two 
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separate  integration  routines  were  established,  one  to  perform  the 
space-angle  integrations  element  by  element  and  another  to  perform 


the  anguLar  integrations  required  by  the  and  operators  appearing 
in  the  integrands  of  the  space-angle  phase-space  integrals. 

An  additional  complication  also  arises  when  the  G and  K opor- 

g u 

ations  are  performed  on  a given  function  of  the  space-angle  variables. 
Although  the  same  function  appears  on  the  right-hand  side  of  Eqs 
and  (22 ) , outside  of  the  integral  as  well  as  within  it,  the  function 
within  the  integral  depends  upon  different  angular  variahles.  As  a 
result,  all  functions  operated  on  by  the  G^  and  operators  had  to 
be  duplicated  so  that  the  equations  could  he  properly  coded  for  com- 
puter use. 

In  order  to  implement  the  G^  and  operations  on  the  computer, 
it  is  also  necessary  to  evaluate  the  scattering  cross  sections 
appearing  in  their  definitions.  From  chapter  II  we  have 
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where  is  the  expansion  coefficient  from  Eq  (19).  In  addition,  we 

have 
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where 


*&,**')*  I (35a) 


ao 

E <5,M(^75L)5(-ii')  (35h) 


Substituting  Eq  (33)  into  Eqs  133a)  and  (3*+a)  yields 


^(t.A-s')  = xT  ‘siX^)(Jxf')[pM  t>  ?(-A-4')l 

/•* 

OO 

^*>11  <r/M(-2ir)[ij(**')-  Pt(  -£*')] 


^(*.**0-  L 'r/‘«(^L)p/(*^') 


(33b) 
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(34b) 


Thus,  by  comparison  of  Eqs  (32a)  and  (33b)  as  well  as  Eqs  (19)  and 
(34b),  we  see  that  « <£s  for  1 even  and  for  1 odd. 

Therefore,  all  the  scattering  cross  sections  can  be  coded  in  terms  of 
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the  known  crojs  sections  which  may  be  either  isotropic  or  aniso- 
tropic in  nature. 

I 

Along  with  the  U and  K operators  and  their  associated  scattering 
g u 

cross  sections,  the  dt- V term  of  Eq  (24)  was  also  coded  as  developed 
in  chapter  II  for  two-dimensional  cylindrical  geometry.  In  contrast, 
the  source  terms  and  as  well  as  the  boundary  term  were  not  eval- 
uated; they  were  merely  set  equal  to  zero  in  their  respective  function 
subroutines.  Several  utility  subroutines,  which  were  written  in  order 
to  further  modularize  the  computer  coding,  brought  the  total  number 
of  subroutines  to  ^5* 

The  final  step  in  applying  the  finite  element  method  to  the  weak 
form  of  the  second-order  even-parity  equation  is  the  assembly  of  the 
global  phase-space  matrix  from  the  element  phase-space  matrices.  This 
procedure  is  complicated  somewhat  by  the  fact  that  the  global  space- 
angle  nodes  as  well  as  the  local  space-angle  nodes  are  located  in  a 
4-dimensional  phase-space.  The  assembly  subroutine  takes  this  into 
account  by  relying  on  two  IE(e,i)  matrices,  one  for  the  angular 
domain  and  one  for  the  spatial  domain,  to  determine  the  global 
space  and  global  angle  nodes.  Once  they  are  calculated,  they  can 
be  combined  to  identify  the  global  phase-space  node  corresponding  to 

r 

a given  local  spatial  node  paired  with  a local  angular  node. 


Once  the  procedure  for  labeling  the  global  phase-space  nodes 
has  been  established,  the  assembly  procedure  can  be  implemented 
as  follows.  First,  form  a null  coefficient  matrix  GA  whose  dimen- 
sions are  (#  space-angle  nodes  X # space-angle  nodes),  or  NSAN  X 
NSAN.  Then  for  each  phase-space  element  e - 1,2...,E  perform  the 
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following  steps.  Form  the  element  coefficient  matrix  |ANSLie 

L Jrnnxmn  , 

where  m = the  number  of  spatial  n^des  in  the  element  and  n = the 
number  of  angular  nodes  in  the  element,  by  evaluating  Gq  (2k)  for 
all  possible  combinations  of  trial  functions,  , and  test  functions, 
. Find  the  global  phase-space  nodes  corresponding  to  the  local  i, 
k spatial- j,  1 angular  node  pairs.  Assemble  the  local  coefficient 
matrix  into  its  spot  in  the  global  coefficient  matrix.  In  abbre- 
viated FORTRAN  notation  the  procedure  can  be  best  summarized  hy  the 
following: 


C 

C LOOP  OVER  ELEMENTS 

C 

DO  5 e = 1 ,E 
C 

C FORM  LOCAL  ANS -MATRIX  FOR  ELEMENT  e 

C 

CALL  GAUSS  (e,.  . . ,ANSF) 

C 

C LOOP  OVER  LOCAL  NODES:  GET  GLOBAL  NODES 

C 

DO  4 i = 1 , m 
I = IE(e,i) 

D03  j = 1,n 
J = IE(e, j) 

DO  2 k = 1 , in 
K = IE(e,k) 

DO  1 1 = 1,n 
L = IE(e,l) 

C 

C ASSEMBLE  LOCAL  INTO  GLOBAL 

C 

GA(I, J,K,L)  = GA(l, J,K,L)  + ANSF(i, j ,k,l) 

1 CONTINUE 

2 CONTINUE 

3 CONTINUE 
k CONTINUE 
5 CONTINUE 
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Of  course,  the  source  terms  would  be  treated  similarly  once  they 
were  defined  in  their  appropriate  subroutines,  "'he  accuracy  of  the 
global  source  vector  and  global  coefficient  matrix  will  ultimately 
depend  upon  the  degree  of  numerical  quadrature  used  in  the  finite 
element  integrations.  When  the  global  coefficient  matrix  and  global 
source  vector  have  been  formed  for  the  assemblage  of  elements,  the 
assembly  process  has  thus  been  completed. 


-tV  Conclusions  and  Recommendations 


The  computer  application  of  the  finite  element  method  to  the 
second  order  even-parity  form  of  the  transport  equation  was  partially 
accomplished  in  this  work  by  the  formation  of  ^5  subroutines  which 
collectively  assemble  the  global  phase-space  matrix  for  solution. 

Because  it  fails  short  of  actually  solving  the  transport  equation,  it 
cannot  be  considered  to  he  a transport  code.  Instead,  the  subroutines 
generated  in  this  work  should  be  thought  of  as  one  of  the  many  approaches 
possible  in  applying  the  method  of  finite  elements.  Of  course,  this 
approach  differs  significantly  from  other  similar  approaches  in  that 
it  applies  the  finite  elements  in  two-dimensional  cylindrical  geometry 
and  allows  the  use  of  anisotropic  scattering.  Both  the  cylindrical 
geometry  and  the  anistropic  scattering  introduce  added  complexity  to 
the  calculations  performed  in  the  subroutines.  Yet,  the  approach  is 
somewhat  inefficient.  This  is  due,  in  part,  to  the  large  number  of 
subroutines  and  to  the  number  of  subroutines  which  had  to  be  dupli- 
cated. Inefficiency  is  also  caused  by  the  repeated  calculation  of 
all  local  shape  functions  for  a given  spatial  or  angular  element 
when  only  one  shape  function  is  needed.  A revision  of  this  procedure 
should  be  considered  in  future  computational  refinements  and  improve- 
ments. 

As  discussed  in  the  last  chapter,  the  four-dimensional  nature  of 
the  space-angle  phase-space  finite  elements  increased  the  complexity 
of  the  digital  computations.  The  doubly  subscripted  trial  and  test 
functions  S.A.  and  S„AT  led  to  four-dimensional  arrays  in  the  local 
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coefficient  matrix  calculations.  These  four-dimensional  arrays  had  to 
be  compacted  to  two  dimensions  in  the  subroutines  because  four-dimen- 
sional arrays  could  not  be  used  in  CDC  TOkTRAN. 

Additionally,  because  of  the  large  number  of  computations  required 

by  the  extra  angular  integrations  in  the  G and  K operations,  the  sub- 

6 u 

routines  in  which  they  are  calculated  may  use  an  excessive  amount  of 
time.  It  is  recommended  that  the  calculations  performed  in  the  operator 
subroutines  as  well  as  the  lower  level  routines  utilized  by  them  be 
streamlined  in  or-der  to  increase  their  efficiency. 

In  conclusion,  although  refinements  and  improvements  in  the  45 
subroutines  of  this  approach  can  be  made,  and  development  of  the 
source  terms  and  the  boundary  conditions  term  accomplished,  a method  for 
the  application  of  phase-space  finite  elements  to  the  anisotropic  even- 
parity  neutron  transport  equation  in  cylindrical  geometry  has  been 
presented. 
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Appendix  A 


Computer  Subroutines 


j 


I 


SUnP('UTINE  OMA 0 ( <M  A ><) 

f; 

4 « 4 4.  « v ««»*  t****.*:'.  .»#•**.  + »««* 

r '■■•trouitit  oiiaj  storf>  <m a y gausstvi  quadrature  rulcs  , * 

n 'I'TRF  IMA  X = THE  HT1PF  > OF  QUAPRATURF  oqtnTF  dFR  QUADRA-  ♦ 
k r T'i-F  RUlr.  < M A X CA  M'lOT  rXCEF0  b . * 

f*  * # »*  + •♦*  ».  4-  V>***-*' 4-**'»**»  ***•»♦* 

c 

rOM-rn/P-'F/I  (c  >,Mr,FG)  ,R(F,50> 

FFA0",  K*AX 

IF  ( K A A X . GT  . 1 ) GO  TO  ? 

DO  1 K - 1 , K 'T  A X 
P c A 0 , 1(0 
T M A X = 1(0 

READ'  , ( (A  (X  ,L  ) ,*1(X,L)  ) ,L  = 1,IMAX) 

1 continue 
GO  TO  7 

7 FPINT*  , ,,M AX  1 M'JM  ALLOWABLE  NUMBER  OF  QUADRATURE  PULES 
JIG  S . ** 

7 ROIJ'  N 
PND 


CUBRC UTI NF  QORULr( T>UL c, JRUL E) 
r, 

rt  > f « t ( 4 4 t a > > -a  4 4 * 4 > » f * 4 » 4444444 

D *“'•  "*  ROUT  J * E PDRULE  DEI  Ep'tl  NCS  WHICH  OF  THF  Av/AHA^LE  ( i>  * 
D ) PUA  Q.-.A  T'JFF  SCTS  TO  RF  USED  FOR  THF  SPATIAL  * 

C INTEGRATIONS  ( JF’OI  NT  > A 10  »'HICH  IS  TO  OF  USED  FOP  THE  * 
D .".'GULAK  INTEGRATIONS  (T’QTNT)  ...  * 

f**«  4 4 4 4 4 4 * * * *•  * 4 » 4 * * a.  4 44»4444**4 

c 

'COMMf  N/ONP/KG  ) , A ( ?>,  F 0)  , 9 (c  , 5 0 ) /T<n/C  ( 50  , D ( SO) 

COPUC  N /THREE/ F (GO)  , c ( « • 0 ) /O'JAOPTR/IPDl  NT,  JPDINT 
FFAO*  , IRUL c» JRU. F 

IF  (If  HLF  . GT  »t . 0°. t JRUL F • GT ,F  ) GO  T3  7 
TPOir'T  = J(IRULF) 

JPCINT  = I ( JnU  LF) 

DO  1 * = 1 * I nGI’IT 

cm  = Ada'MLF.K) 

0(0  = F)  (IRMLF.K) 

1 CONTI NUF 

DO  ? L = l, JPDINT 
P(L)  = At  JOULSU 
F(L)  = °( JRULE|L)  ' 

? CONTINUE 
GO  TO  <* 

o PRINT*  ."MAXIMUM  NUMBER  OF  QUADRATURE  RULES  IS  b 
U RETURN 
FNO 
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* ■ • /.r 


( 


FUNCTION  0(1) 

C 

4»*-»4*«»*  > *•****  ' »«  4V»44»444 

r ” ’ I M c T 1 P N!  0 CALOUlAYF'  ri)r  PUAOfATURf  POINTS  FO'  THF  * 

0 T^nif.'T  ^l.'A ORA  .TURF  RU.E  * 

r»  *♦****«*<►»  #44*4*  »»  ***»v*»  4*4 

C 

COPRFN/TUO/C  (F  0)  , o<  • 0) 

° = r (I) 
rEUiF  t; 
cNn 

F'JhCT  ICN  W(I) 

r 

C r I NCTON  W CALCULATE?;  TUC  CUA  ORA  TURE  HEIGHTS  FOR  THE  * 

0 IP  01  NT  Ol'AOP.nTUFE  RULE.  4 

0 * ***>f44»'r»*«.44444x»44*44444»» 

C 

CORMCN/TW0/C(.;O)  ,0(30) 
w = cm 
P r T U r N 
END 


F'JNCT ICN  00(1) 

C 

C*  *..*4**.»4»  ».  44  * *4r»4«  4*4*444 

C r'INCTIO!'l  C A L CUL  A Tr  S THE  OUADCATU^  DOIMTS  FOP  THE  * 

C JPOINT  CP  A ~jT.  T'J  F E RULE.  * 

C*  4444*4  i-444444444*  * 44*4444444 

C 1' 

C0FMCN/TH9EF/E  (30)  ,c(Fn) 

CO  = F (I) 

RETUF N 

fnc 

■* 

FUNCTION  WW(I) 

r 

p«  4 4 4 4 4 4 •. 

C r-inCTICM  WVf  CALCULATFS  THF  QUAORATUvE  HEIGHTS  FOR  THF  * 
r.  POINT  OI’ADRAT'JP E RULE.  * 

0H  4 4 4 * » 4 <-4  4 V 444444^4  444  4 4 * » 44* 

C 

COf'UCN/THRFr/C(aC!)  ,r(*;o) 

WW  = F ( I ) 

FE7UCN 
ENO 


38 


c'U<ixt  UT  I ur  AES  M^L  (M'WP,Max‘t»  MAXHU.  MAXCHI  , NSN,NSAN,N,M, 
m.,CR,C7,  CArVJ,CCHI  r E,  JF,SF,  AE,  G\,  GS,AS~P,At4jr-) 


c 

c 

c 

c 

r 

C 

C 

C 

f 

C 

C 

r 

c 

c 

c 

r 

c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

c 


*■»-**.-  v***!*1'  « t « < r i « > f i * » # * + 

SlWOUntF  ASSVRL  r»rT'-'°HINrS  TJ E "STTrrNESS"  MATRIX  AnO  * 
•■l.OA'V*  'If  C’CK  rnx  <-;*H  GPf  rE-ANGLE  DHVSE-SP«CE  ELEMENT  * 
AMO  A SG^PlES  nU  T H INT0  furi^  glo°AL  C 0 J N T Ef  =>  A R TS  . * 

A^GM^L  lr  C A l L r O N IT  H TUr  *-AXIMUM  NUMBER  OF  R,7,MU,  ANO  * 


THE  VALUES* 
OF  SPACE-  ♦ 
’SHAPE"  F UNC * 


KL 


rwi  ^ E C.H  Li'JFG  AMO  T h ■ I P CprkEGpoNOTWG  ARRAYS, 

Fir,  THE  tU-P-R  C;F  A I “ U L *V 0 f'OOrS  ANO  I HE  NUMBER 
AMGl r PHASE-fPACt  M0:rG,  fur  MJHOER  3F  LOCAL 
ATOMS,  At"!  Tt-E  IE  E J-  GLO^Al  NOOF  MATRICES. 
la  gel f(i  rr»-Mou  ol  j r.y  htlow  is  used  to  transmit  the 
lap  “'-ch  a r pa y r io  c.uigtt o"  gausg?. 

L s °El  - C f Ot'HOtl  " L 0 C'<  EL-hf»;T  PASSES  ANGUL  AR  ELEMENT 
FUNCTION  gauss?.  * 

t.'CMEt  CL  ATIJ PR  F OR  A ® >UMrNTS.  . . . * 

H ■*  XP , ha  XT  » H v M’ ) , M f X Oh  7 - MAXTMJM  11  0P  MFSM  LINES* 
NaN  = t'UMnR5  3-  NROFS  IN  THF  ANGULAR  DOMAIN 
NS  AM  = NiJHOFR  OF  VPACE-aNGLF  PHASE-SPACE  NODES 
m = MUM  REP  3r  SPATIAL  TENT  F JNC T I ONS ( L = L INF A P , 
6=0U\0-;f.TIG,  E“C  . ) . 

N = NUHP£F  pr  ANGULAR  TENT  FUNCTIONS 
IF  = GLOOA.l  SPATIAL  node  matrix 
JE  - GLOBAL  ANGULAR  NOn  E MATRIX 
SE  = c P A T J V L rL  E"rNT  NUMOfR 
AE  = ANGOLA  * CL  FH rNT  NUMOc"‘ 

GA  = GL  01  A_  ’*  G T T rr  N ESS"  MAT?IX 
GS  = GL 00  AL  " L 3 A 0"  VEGTQP 


A NGU-  * 

* 

TO* 


COt.uc  N / EL  E MFNr  /< L 

COMMCN/HILOW/ALO  (3  0)  , AMI  (50)  ,CMIL3<5  0 > fCHIMi  («,o> 

Pit- FUSION  CR  (*1AXC  ) ,G7  (UAX7)  , CAMU(MAXHU)  ,CCHI  (MAXCHl)  , 
$ AuSF(mn,mn>  ,A  YSG(  MM)  , 1 £ < L> 0 ,1B*  , JE  (50,  la)  ,GS(NCAN>  , 

5 GA  (NS AN  ,NS4N) 

INTEGER  3E  , A E 


C 

C TNTTIAL17E  GLOnAL  MATRICES. 

C 

PO  4 K1FST  = l ,M 

•a  = IE  (SEXIEST) 

00  7 LTEST=1,N 

L = JE(AF.LTEST) 

LXN  = LTEST  ♦ (KTFSf-l ) *N 
XL  = L ♦ ( K- 1 ) *N c N 
OS (Kt  ) = 0 . 
no  ?.  ITRJAL  = L,m 

1 = 1ECSE,1T»IAL> 

PO  1 JTK’I  AL=  1*  N 

J = JE  (Ac,  JTRI  AL) 

JIN  = JT->IAL  (ITRIAL-1>*M 
IJ  = J ♦ <I-1)*MM 
GA(<l,IJ)  = 0. 

1 CONTINUE 
? COf.T  I M'E 

7 CONTINUE  ,Q 


r»  o r»  .-»  o o 


A CONTINUE 

xaxrni  = m ax  r - : 

MAX7E1  = M ^ V 7 - l 
MAXMI'mJ  = "IX  HI  - 1 
*JAXC»'M3  - "AXCUT  - t 
mAXSF  = f’AX7“i* 

H A X A { =■  m ft  y^iM  l*  “ A VOHM 1 

RAIL  HlLO<  ALO, AH* , OHJ  L 0, CHIHI , MAXI J,  M AXCHI,  MAXMUM1 , 

J MI  yCHMl  ,»AX  AR^AHl^rcHl) 

LOOP  OVH  EACH  t L E HE * T IM  RHASE-SPACE. 

PO  1?  I = 1 » HA  X 7 Ml 
7LO  - C7 ( 1 ) 

X H I r 07(1*1) 

DO  11  JM.-AV^Ml 
PL (i  - CP(J) 

PHI  = C<(J*1) 

PC  1L  *=1 , wAXlHMi 
pp  ° l = i,  m Ar'ijMi 

*<*L  = L •*■  ( K—  1 ) ■*  H A X M' JM  1 

CALCULATE  THE  ELrH~'|r  AL  PHASE-SPACE  "STIFFNESS"  MATRIX. 
PAIL  GAUSS  <N  ,M  ,'1 v'  f Re  0 i r H I , ?l  0 , 7'*T  , A L 0 < Kl  > , A HI  (XL  ) , 

t r hi hi (xl) ,:hilo(xl) ,ansf,anss) 

c 

C LCPP  L’VE-*  LOCAL  \OOrS,  GET  GLUPAl  NODES. 

c 

SE  = J ♦ « I- 1 ) X7H1 
AE  = XL 

PO  A XTEST=1 ,N 
x<  = I E (S f »K  TEST ) 

PO  7 LTEST=1»N 
LL  = JF(AE,LTEST> 

L KN  = LTEST  ♦ (KTRST-l  )-*N 
XKLL  = LL  ♦ CKK-1  > *N«N 
F-SCKXLL)  = GS(XX'.L)  ♦ ANSG(LKN) 

DO  6 ITRIAL=1,M 

II  = If(SE,ITRIAL) 

DO  b JTPIAL=1,N 
JJ  = JE( AE.JT’IS.) 

JIN  = JTRIAL  «•  (1T°IAL-1)*'( 

IIJJ  = J J ♦ <II-1>  * JAN 
C 

C ASSE*‘QLE  LOCAL  INTO  SL 03AL . 

C 

GA (Kxll  » T I JJ)  = GA (KKLL , IIJJ)  ♦ ANSF ( LKN*  JIN) 

? CONTINUE 

t CONTINUE 
7 CONTINUE 
A CONTINUE 

9  CONTINUE 

10  COnTtnUE 

11  CONTINUE 
1?  CONTINUE 

RETURN  Ln 


I 


suokc  in  i ne  HiLoo.Lo,AH',rLOfCHii'r,,nymi,‘«£xCHi,*rx*ijMi, 

S M/yrHMi,MAX.0*.MU,3C»‘I> 

c 

0 * * 4 ' 1 4 4 * ♦ * *«  „ »«*♦**.  • t » * * * * * 

C St'OROUTIt  F HILO  CU"J.«Tr«;  THE  A L 0 ♦ M T * 3 . 0 ,A  NT  CHIMI  VcC-* 

C TORS  FOR  f/.Ch  FI  F'l  F H I *4  Tr<r  ANGULAR  'OHIN,  • 

- » t « t • 

c 

OIMt'f.fTON  ®L  0( M®X)  , AJI  (MAX)  , OLOMA-X)  , CW1  HI  (MAX)  , 

S Of  vic-'A  ymi|)  ,"mt  (‘'AVCHI) 

C 

C LOOP  OVFi  fc\PH  A f,3  Ul  - 5 ILC>'ENT  ftuD  SOU  T ^F  H 1 AMI  LO 
C VECTOR. 

C 

00  2 I=1,MAXCHM1 

no  i J--1 , MAXMIJ’11 

1ELEMNT  - J + (I -1 ) * *1 A VMf IM1 

A L C ( 1 ELEUNT)  = 0 1 M t)  f )) 

AH1  ( IELE'-*MT)  = CLMIK.I  + I) 

CLOUFLc,UIT)  r CrHIC  ) 

••  OHiwj  (JELEM'IT)  = C0;-»I  (Tfl) 

1 CONTINUE 

2 CONTINUE 

k El  Ur  N 
END 


SUBROUTINE  GAUSS  ( M , M , MU , RLO , RHI , 71 0 , Z HI , A L 0 , AHI , OHIHI , 
JCHILC, ANSF, ANSG) 

C 

CR  ♦ ■*■***  **♦*  r ****♦♦  v » ♦ 

C INTEGRATION  of  F(R,7,  MLI,CH])  .OR.OZ.OMU.OCHI  BY  GAUSS  IAN  ♦ 
C QUADRATURE.  * 

C*  • *4«+4*»44>-»4444****4*f4**».4* 

c 

REAL  1NTGRQ1 

COMMON /SPAC  F/XI<tFTU  /TESTS/KTEST  ,LTEST 

COMMON/ GUAOPTS/I3OINT,  JPOINT /T PT ® L 3 / 1 TRI AL  , JTRIAL 

DIMENSION  4NSF(MN,MN>  , ANSG(MN) 

C T‘IITIALI7E  ANSF  AN  Q AN  S3  \/A  RIAnL  ES . 

C 

DO  1 KTEST  = 1,M 

CO  1 LTEST  = 1, N 

LKN  - L T E 3 T ♦ (KTEST-1)*N 

ANSG (L  KN)  = 0. 

DO  1 ITRIAL  = 1 »’1 

DO  1 JTRI AL  = 1 , N 

JIN  = JTRIAL  «■  (1TRI AL-1) *N 

ANSF(LKN, JIN)  = 0. 

1 CONTINUE 

C 

C SET  UP  COORDINATE  MAPPING  OF  ELEMENT  IF,  I = 


o o o o o 


CALL  Mf  M^LO,  RUI*SFA3  ,ZL0,ZHI,7FAC,  AL  O.AHI » AFAC.CHILO, 

t r.*UM4,CHIF 

FAC  =■  T F->  C ‘ ■’"AC  * AFAC  * CHI  FAC 

r 

r FVALI  A T r NEW  VaR  T".lLrS  WHICH  A Rc  EX^RFSSEO  IN  TFRHS  OF 

C THE  l EGEin  lr  ROOTS  . 

C 

DO  3 Irl1I'»OI,<T 
XII  = 0 < T > 

TO  3 JsliIPOINT 
FTfJ  = O(J) 
no  T K=1 * JPOINT 
XI  K = GO  ( K ) 

HO  3 L = 1 » J^OINT 
FTAL  = OO(L) 

WEIGHT  = W(I>*U(J> *WW(K)*WWCL) 

C 

r SET  I’P  TER  IS  TO  3E  INTrGRATFD  A*>n  CARRY  OUT  INTEGRATION 

r,  «Y  CfLCULATlNG  A NS  = SUM  OF  W*F(XI,FT4)  . 

C 

CO  2 KTFST=1,M 
00  ? LTEbT=l,N 
L*N  = LTEST  ♦ (X 1 EST  ■*  1 ) * N 

ANSG(LKN)  = ANSS(.KN)  ♦ WEI  GHT  * TNT  3RD  2 (X I K , ET  AL  ,X  II  t 
5 FT  f J) * F AC 
DO  2 ITRI4L  = 1,M 
00  2 JTR1AL  =1,N 
JIN  : jmi A L (TTRIAL-1)*N 

ANSF(LKN, JiN)  = ANSFtLKU, JIN)  ♦ WF I SHTfc I N T GRD 1 ( X I K , 

3 ETAL.XII ,ETAJ)  *FAC 

2 CONTINUE 

3 CONTINUE 
RE TUP  N 
END 


SUPROUTINE  MAP (WLO,WHI,WFAC,XLO,XHI,XFAC,YLO, YHI,YFAC, 
J ZLD/  7~tl  , 7FAC) 

C 

C*  ' 

THTS  SURF OUTINE  COM°UTFS  THf  SCALING  r ACT ORS  WFAC,XFAC,  • 
Yr  AC » ANC  7FAC  USED  TO  MAP  COORDINATES  IN  THF  W,X  L Y,7  * 

“LANES  If'TC  THE  W',X»  * Y*,Z»  PL  ANTS.  • 


COMMrN/MAPSPAC/A,OfC,n/MAPANGL/F,F,G,H 


A 

= 

(WLO 

♦ 

w-m  /?. 

fl 

= 

( WHI 

- 

WLO)/?. 

C 

= 

(XHI 

XL  O) /? . 

0 

= 

( XHI 

- 

XLO) /?. 

F 

= 

( YHI 

YLO)/?. 

F 

r 

< YHI 

- 

YL0)/2. 

G 

s 

< 7 HI 

♦ 

ZLO)  /?. 

H 

s 

(7HI 

- 

7L 0) /2. 
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WF  AC  = B 
XFAC  = 0 
YF  AC  = F 
7FAC  = H 
f<F TUf  N 
FNO 


-•FPL  FUNCTION  INTGRD1(YIK,ftal,XI:,FTAJ) 

^ l ***.*.*>  • ♦ ♦ * 9 ♦ * * * * ¥ 

0 SUBRr'JTIt  - IUG'.Dl  Of  LC  ILf  CS  THE  1 I T F “>  R A N 0 f t-  THF  j T?  FAM* 

0 T -w,  f.-  «.  CA’TF  UiH  3HA  ~p  sr^CE  RAT!-V  Th  T^F  UEAK  cO*m  0c* 

0 "r\IE\  “-Pa  U7  Y SFLON)  Of  TER  TRANOD0RT  EQUATION.  * 

C '••nKOUTIt  E IUTTAOl  TH-  FOLLOWIJG  fjnction  suo/r*s  * 

0 P.KIJ  OOCLTST  BNDCCS  * 

0 r,  C-  ODFLTFL  OPELTR2  * 

0 TRIAL  TEST  TF1AI?  * 

rj  r a*****^*  »»»»«»***» 

p 

CX1  Fi  f,AL  3°l  F LT  L,0^cLTB?tTrIAL,  TRIAL? 

Fl  = FK'MOJElT  •’.L,0"'FLTR?,XlK,ETAL,XiI  , ETA  J)*OOELTST 

t (yi«,:tal,x it, eta j> 

C2  = r,r.(TRirL,TRTAL2»XIK»ETAL,Xir,ETA  J ) » TEST ( X I < , ET AL , 
> Xil,ElAJ) 

f-»  = i&ial<x:*:,ftal,yit,etaj)*T':?t<xik,et4l,xii,etaj> 
i *■-  r.TFCi ( xik.fj  AL,xT:,Eirj) 

INTGf  Pi  = Pi  ♦ p?  f FJ 

r-FTURN 

FNO 


FUNCT  ION  r<<U(FFT,cCT?,  XIK,ETAL,XII  ,ET  AJ) 
r; 

I',*  ***4*-**-**»  * « ♦»**•«•  I k********* 

F FUNCTION  Ri'l.l  OETFRHINES  the  VALUE  O'7  THE  OPERATOR  RU. 

PV-  ♦ * + *•»  »*♦■*♦♦“■**¥*  k ♦ * v * 4 «► 


^YTERuf  L Sir,ilA<UfFCr2 
= p (XIK) 

7L  = 7(PrAu) 

A»UI  = A»‘IJ(XJI) 

PH  1 J = CHKFi'AJ) 

FKl'  - FCMX1K,  ETAL  ,XII  ,ETAJ) /XT(P<,7L)  ♦ 
G/  USu?  (SIGNA<U,FCT?,RK,7L,AH'!T,CHI  j) 

ff'.'FU 

FNp 


k 
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C'C  ' I Fur. 'T ICN  INTOF  0?  ( y [ < , T AL  * yI  I , ET  A J) 

n 

Q*  * * ^ ♦ -f  * 4 ♦ ♦ I * * ♦ ♦ * * + 

r O'lORCU'  II  c IlTJfT?  «V  t C'lL  ATS  THc  T'i[rr)<ANOr  QT  THr  * 

o <?'»iircc  vrno‘  rr  t te  we\k  rj  p n of  tip  ,,gven*,-3apity  second 

f'  nROFt  1W  |.f  nr  r Ff>UATIO‘J,  ♦ 

C 'L1  •■’pr’J  | II  F Ii.l.PL?  JTS  r N r FOLLOW!’''-  FJNCT10k,  su"»/r»s  * 

C SG  TEST  * 

0 RKU  OPELTS*  * 

0 011  * 
' F')NCTTON  SU  ’.'HEN  OPCMT:"'  uN  nv  F'Jf '0  T 1 0 ‘J  RKU  WILL  REQUIRE* 
r ,,nOIrICA*  I ON  ( > PE  f J*  r TI3NC  0 D F L T r ? AMO  TRIAL?)  • * 

r>v  ♦«**♦*<**  *»«*«»*#**• 


r x’  “r  K L TJ 

rl  - I G(x  t<  , rT  ul  ,XIt , PT  J)  A TpST(  XTK,  FJAL,  XII, ETA  J) 

F?  = R>''U<F'l,Xi  SrTAL,XII,ETA  J) ‘OTTST  (XI<,ETAL,XII, 

* FT  A J) 

7NTGT  0?  = c\  ♦ fo 
F ET  uf-  N 

FN[.  *. 


C,J^CT  ION  GG  (FGT  , FCT?,  VK,  ETAL.XII , ETA  J) 

r. 

i - « ♦/•  + »»♦*♦»*»##«***.  * * # * » * * V * ♦ * 

0 FUNCTION  GG  ?FTFr  •*I,|F«5  THE  VALUE  Or  THE  0 PEE  J TOR  GG.  * 

*«.*  «.,  + »*  *■  1 » « 4 * + * ^ * *■»**■*¥■**  + * 

FXTFENAL  Si  G'-lAG  , PC  T ? 

RK  = E(XK) 

7L  = 7 (ETAL  ) 

A’lUT  = A-'J(XH) 

E'JIJ  = CWI  ( ET  A J) 

GG  = XT  K » "'L ) F?T  (YIX,  ftal,  XII,  ETAJ)  - G MJSS  2 ( E IGMAG  , 
2 Ff  T ?,E  Yj-'L  , A’HII  ,GHI  )) 

PETUPN 

fn& 


FIJNC*  ION  XT(  R,?) 
r 

0*  ♦ ********»..  **444.,  ,,4*444444,. 

r FJNCTICN  XT  CALCULATES  Mr  ANISOTROPIC  TOTAL  CROSS  SFCT.  * 

r 

XT  = 1. 
fno 


* >'V 


FUNCT  I(M  OQE  L r (XI<,r-AL,XII,ET\j) 


c 

r*-«**'-'*-»'**»*  + '‘*«***'i'*»**»*4-»***. 

C FUNCTION  OuELTRL  :A!.:JlATFc  The  omega  tdi  OFl  trial  term  * 
C (THE  ST  vr  A Ml  MG  TRH)  IN  CYLINDRICAL  COORDINA  TE3  • *■ 

£4*  *44  4444444  * 4 4 4 44  . 44  444444444 

c 

co>'M'N/Hcm"s/M  ,Nj/TTmr/i,j 
T3  = TRILL (<i<  « “ T A L > X I t t ETAJ) 

call  rrRrv(M,M,i»J»r>i  ,n2*nT) 

RK  = ‘ (X1K) 

AH'JI  = A *•'  U < x 1 1 ) 

CHI  J = CHI  ( FT  A J) 

A - 1.  - A C II*  vo 

9 - rO.'T(A)  / »■< 

OOFLTRL  = r>‘C3;CHI  I)  * ( (R<»01)  *T?)  - R* c 2 * AM'JI*D3 
P E 1 1) c N 

FNn 


FUhCT ION  OOELTP?  <XIK, r-AL,VII,FTA J) 

r 

C*  * ♦ * * 444444*  »*  + ***-•  » *•  ♦*  4444444 

C FUNCTION  0CFLTR2  CALCULATE  OMEGA  TOT  OFL  TRIAL  TERM  * 

C <THF  STEMMING  T CRM)  IN  C VI  INDFICAL  COORDINATES  FOR  USE  * 
C IN  FUNCTION  GAUSS?.  * 

0 * * 444444*4*  » 4 * 4 4 4 4-1*  * 4 * 444  « 444 

C 

COMMf  N/NODcS /M »N/TRIAL~/I , J/LNOOF/  J J 
T°2  r TRIALS  (X  K, RIAL  ,VII»ETA  J> 

CALL  0ERI7?(M,N,I  , J 1,03 ,02,03) 

PK  = P ( X T K 1 
AMUI  r AMU  (XII) 

CHIJ  = CHI  (ETA  n 
A = 1 . - AMUT**? 

0 = f CRT  ( A ) / R’< 

0 0 C L T K 2 = n*C0S( CHI  J) * ( (RK*01)  ♦ T R? ) - 9*?2  ♦ AMU1*03 

rsTurw 

FNO 


< 


FUMC • ION  OOELTSTfXIK,  FTAL,/xTf pT;  j, 

C 

C*  * 4 4*444.  *4**44444  ->  I*  ♦»  4*  4*4*4 

C FUNCTION  rnCLTf.T  CALCULATE0  THF  omega  tot  OFL  TE°T  TERM  * 
C (THE  STRTAmims  T°RM)  In  CYl  INDrIGAL  COORDINATES.  * 

£ 4 * 44*4444  %4*4  4444  4.  * 4 4 * 4*.  44**4 

c 
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I 


r ' ■ 


c r u o i'  o o o 


rr;  ■>.,  / I , J/r'’ST<  /<  ,L 

T ,;n  - i C jj  (XT<  r \L,  *1 7 , ETA  J) 

rm  nrxiV(M,M,<tL  , 11  ,r2,n?) 

<.'W  = ^ < X T < ) 

A*"  MI  = A-MI(XII) 
ruiJ  = CHI  ( F T A J) 

A = 1.  - A “I ) I » ^ ? 
r>  = c or  T ( A ) / ( < 

rnaTCT  - r*cds(  :hi  j)  * r ( *<■  on  »■  tsd  - r*o?  + amui*d.t 

«-r  T'i»-  n 

’NO 


FUNCTION  T»I1L  (XT*, FT  A!.  ,XII,  ETA  I) 

c 

0*  •♦♦**'***'»»'*'***»»-*f*.^**»»***4 

n F (NOTION  TRIAL  CALCULATES  THE  VALUE  D-  THE  TRIAL  SOLUTION* 
C ; ( T ) * A ( J)  . * 

Q*  • •♦■♦•♦♦♦♦♦•♦*****^»*4**4*¥4*# 

0 

COMMON/ TT  IALS/I,  J 

rOMMC'N/TfNT  1/  S,  01 0X11  , 0T0ETA1/Tr  Nr?/A,QTDX12»DT0ETA2 
DIMENSION  S(16),A(16),nTDXIl(10),OTOXI2(16), 

S OTDETAl < IB) , "TOETA 2( IB) 

CALL  TFNTENClXIKjETALtStOTOXIljOrDErAl) 

CALL  TENTFNCIXII,  ETAU,  A , OTCX 1 2 , DT 0 FT A 2) 

TRIAL  = S ( I ) * A ( J ) 

RETURN 

END 


SUBROUTINE  SCALE ( N,DTDXl , OTOFTA, AMJ- AC , T HI E AC , DT DC HI , 
5 OTDMU,T,CHI  ,XSTN) 

C 

0 * r ♦ »*♦♦*«*♦*****♦**»********♦ 

THIS  SUBROUTINE  CALCULATES,  FOR  A GIVFN  ELEMENT  E,THE 
DERIVATIVE  Or  THE  "TEST  FUNCTIONS",  MI),  WITH  RESPECT  TO 
CHI  AND  NU  BT  SC A L INC  THE  VALUFS  OTOXI  ANO  PTOETA  BY 
(DFTA/DCHI)  = 1. /CHIRAC  ANU  (DXI/OMU)  = l./AM’JFAC.  AD- 
o: T IONALL  Y , XSIN  = THE  ? A°TIAL  DERIVATIVE  OE  f « I ) * 

SIN  (CHI ) WITH  RESPECT  TQ  CHI  IS  COMPUTED. 

AMUEAC  = SCALING  -ACTDP  SUPCL IPD  3/  SUBROUTINE  MAP. 
CHIRAC  = SCALING  FACTOR  SUPPLIED  3Y  SUBROUTINE  MAP. 

0*  • ♦ ♦*•♦«*♦******♦**♦»*♦*►***♦* 

C 

DIMENSION  DTOXI <N> , DTD ETA (N> , OTDMU ( N ) , OT0SHI (N) , T (N) , 
l XSIN(N) 

"0  1 1*1, N 

nTPCHI  (I) =3TDETA (I ) Ml  . /CHI FAC) 

OTDMl’  ( I ) = DTDX I (I)  * ( 1 ./AMUFAC) 

XSIN( I) -DTOCHI ( I) *SIN (CHI) *T (I)*:D5 (CHI) 

1 CONTINUE 
RFTUF N 
ENO 


\ 
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u C'  o o o c~  o o c: 


SUB PC  UTINE  TENTFNC  (XT  , ET  A , T , 0TO< I , DT  OET  A ) 

DIMENSION  T (t  ) ,010X1  (A  ) , OTDETA  (4) 

C 

C*  • 4444***44,  4*4444*»44**44..44* 

C MTS  SUOf-OUTINr  CALCJLATFS  THE  VALUES  C0R  THE  TENT  FUNC-  * 
C T T 0 NS  AND  THE  I,<  D EKI  V A T 1 V Fr  ON  THE  OASTS  OF  A CANONICAL  * 
C ELEMENT  WHICH  IS  A LINEAR  SFRENOIPIT  i RECTA N3LF  WITH  * 

C N"DES  AT  (XI,FTA)  = ( - 1 , - 1 ) , ( 1 , - 1 ) , ( 1 , 1 ) , AND  (-1,1).  * 

PI*  4***44*44»*44444*,  4*4,  ■,*,*** 


CALCULATE  THE  TENT  FJNCTIONS  T(I). 

T ( 1)  = ( l.-XI )* ( l.-ET A)  /*, . 

T ( 2)  = < 1 .+  XI  )♦  ( 1 . “F  r A ) /A. 

T( 3) = ( l.+XI) ♦( l.+ETA) /A. 

,T(4)  = (l.-XI)*  (l.fETA) /*♦. 

"ALCULATF  THE  DERIVATIVES  OF  THE  TCI)  WITH  RESPECT  TO  XI. 

PTnxi <i)=(FTA-l.) /h. 

DTDXJ ( 2)  = -DTD  XI ( 1 ) 

OTDXI <3)=(l.*ETA)/4. 

DTDXI  (A  )=-DTDXI ( 3) 

CALCULATE  THE  DERIVATIVES  OF  THE  T (I)  WITH  RESPECT  TO  ETA. 

OTOFT A (1  ) = C XI - 1 . ) /4. 

OTDETA  (2)  =(  -l.-XI)  /4. 

OTCETA (3) =-OTDETt (2) 

PTOETA  (4 ) = -DTD ETA ( 1) 

RETUPN 

ENO 


FUNCTION  TEST(XI<,ETAL,XII,ETAJ> 

C 

C*  * 4444444444*444444,44*4444*44 

C FUNCTION  TFST  CALCULATES  THE  VALUE  Q-  THE  TFST  SOLUTION  * 

C S(K)*A(L)  . * 

C,*  * 44444*4444444444*,  44**44444  4' 

C 

C0MMCN/TF3TS/K,L 

COMMON/Tf NT  IX S, DTD XII, PTOETA 1/TEVT 2/A , OTDXI 2, 0TDETA2 
DIMENSION  S ( 16) ,» (16)  ,DTOXI 1(16), DTD  XI 2(15) , 

$ DTDETA1(1»S)  ,DT0ETA?(16) 

CALL  TEN!  END (XIK, ETAL , S, OTPXI1 , DTD  FT A 1) 

CALL  TFNTFNC(XII, ETAJ, A , OTDX 1 2 , DT D FT A 2) 

TEST  = S<  K) •A(L> 

PFTURN 

ENO 
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SUSEC  UT1\E  0ERlV<‘Vl,I,J.01,02f"m 

0 

0*  ‘ «*».»**.***»****¥v'-f»*«  »**»■** 

C Cl  CALCULATES  THE  VUUF  PA"  T I AL  ( R»  FZ  T ) /F  A RT  X AL  ( R ) 

C 4 ° E A R I L ' r I N THE  3 1 R 7 A M I ‘ | G TtRH. 

p 0?  C t LrUL  4 T ES  THE  V AL  Ur  P A- T I A L < f 3 T»  f,  I N ( 3 HI ) ) / ° A PT I AL 

0 ~jT  APPEARING  IN  THE  3 T >r A “ING  TrRM. 

C 0*  0/  L CU_  A T r3  f Hr  V 4 L Ur  PART  IAL  ( r3T)  / PAR  T I AL  ( 7,)  AP- 

C 3"  A RING  I N THE  STPCAMM3  T rRM . 

C * ♦ ♦♦+*♦***.*»•♦**#  '.»*»*  + <«» 

c 

CO*HC  n/tfni  i/s,o:ox i,  oeoeta/ten~>m,  paoxt,  daoeta 

CCMON/’-’f  P3°A0/Ar  , i\rrAC»Cf7FAC/MAc6H3L/c»  A- AC,G,  CFAC 
nlMFt  r.zor  ISOXi  < 1 L ) , l^ETAdb)  ,03D7<  ib)  ,0300.(16)  , 
j Pf  npl , ( 16)  , DA  ■'XT  ( 16)  , OADETA  (!  ?)  ,DAOCHI  (1  6)  ,S  ( 16)  , 

3 A(  !ii)  , 3 A 3 HU  ( 1: ) , )A0CIN<16> 

CALL  SCALE (M, OSOXI , T30 pT A , PF AC , 7 7 A 3 , 0 S07 , 3 SDR , S , 7 , 

i crosiu) 

CALL  SCA_E(  l!,3AD*[  , DAD  rT  A ,ArAC,rr40f  0 A DCH I , OADHtJ , 

E A . CHI , DA  f)  SIM) 

D1  = A ( J>  ♦OSOR  (I) 

P 2 = S ( I) * 0 AOSIN ( J ) 

03  = A ( J)  * US07  < I ) 

PFTIJFN 

FNO 


SU9F  CUT  INE  Or.  RI\/?(N,'J,I,J,01,02,03) 

C 

*■*♦*♦*♦»**  ***•*>  ,»**  ***.**** 

C T-dS  SUP»*  CUT  I ME  IS  NSFO  OY  FUNCTION  DO  rL  T R2 . * 

r nl  CA  Lr U. A T ES  TUr  VUUr  P ART  I AL ( R* FQ T ) /P A R T i A L ( R ) * 

C APPEARING  IN  THF  STRF4KEN6  TERM.  * 

C H?  CALCULATES  T H r VALUE  P APT  J AL(  F C T * S IN (CHI) ) /PARTIAL* 

C CHI  APPEARING  in  the  'TEAMING  TFRH.  * 

r 03  C t LCIILAT  rS  T Hr  V A L LI  F PART  IAL  (FC  T ) / PART  I A l.  ( 7)  « P-  * 

C FEARING  IN  TM£  STREAMING  T^RH.  * 

C “ #■**•»■*'**  + **»*•**■**•*♦ 

C 

COvmc  N/T ENT  1/S, OSOXI,  Q^OETA  / TEN’’ T W.D  / A , OA  OX l , OADETA 
CO.HMf  N/MAPSPAC/AC  ♦?- AC,C,ZFAC/MAPAN32/X,  APAC,Y,CFAC 
01  MR  SI  Ci.’  JSOXI  ( 1 if ) , OSnETA(lG)  ,030  7 ( 10)  ,0SDR(16>  , 


i OSOSIN(IF)  , 04 ' X I ( 1 G ' , OAOE  TA ( 1 j • .OAOCHI ( 16)  »S(  16)  , 
$ A ( HO  , OAOHU  (1  s ) , OAPSINda  ) 

CALL  gcalEL  M,DSO/I , 13PrTAfRFAC»7rAC»DSOT,OSOR»3fZ» 
i DECdt.) 

CALL  S C A _ E ( N , OA  OX  1 , 0 A 0 CT  A , f F AC  , 0-  A C , 0 A np  H I , HA  DM’J  , 
i A , CWI , OADSIN) 

01  = A ( J)  "OSORd) 

P?  = S(I) ♦OAOSIN ( J) 

.03  = A ( J) * 0 S07 (I) 

REI Uf N 
FNO 
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5 

J 


FUNCl  inn  TRTAL2(XIK,ETAL»XII,ETAJ) 

C 

p*  •»♦♦»•■»«*»*  t*****...,**,***** 

r.  r'l  NCTIGN  TRIAL  2 CALCULATE^  T Hc  VALUE  OF  THE  TRI A L SOLU-  * 
r T ' ON  > ( I ) 4 A ( L QC A L)  T?  OE  NUMERICALLY  INTEGRATED  IN  * 

r FUNCTION  GaUSS?.  * 

44*-.4*4***'--*44*4  »,  4*.  44  V 44 

COMVC  N/T..IALS/  I,  J7L  nOE/LOCAL 

COMM(  N/TrN|Tl/S,9T  DXI 1 ,9T0ETA  l/T=gTTRO/A,DTPXI  2,  Q THETA  2 
PIKFMilUN  3 (It) ,f (IS) irTDXIl (lfc),0TDXI2(13) , 

S OTOE  TA 1 ( 13)  , nTOETA  2(  l»'.l 

CALL  TEN' FNC (XI<, ET\L , DTOXIl f DT DET A 1) 

CALL  T rNT  FUC(YI£  , E7\J,A  , OTD X 1 2 , rm E T A 2) 

TRIAL?  = 3( I) *A(L0CAL> 

RETUT  N 
END 


FUNCTION  SIGMAKU(R,7, AMUI,CHIJ,AMUK,CHIL) 

0 

C*  ' *44444444»  444*44  *4  44*  44444*4 

C THIS  FUNCTION  EVALUATES  T N c ”Kl»"  MACROSCOPIC  CROSS  SECT.  * 

C*  4444*4444.  4444*44444*4*4*444 

c 

AMUNCT  = UZERO(A''UI,CHIJ,AMUK,CHrL) 

PI  = 3 • 14139  26533 
SIGMA  KIJ  = 0 . 

L = L FOPrL< IOUMMY) 

IF  (L  . FQ.O)  RETURN 
00  1 1=1, L, 2 

SIGMAKU  = SIGMAKJ  ♦ (?"  IM)*XS(k,7,I)*°OLY  (I, A MU  NOT) 

$ /t./PI/(XT (R,7) -XS(P,Z,I) ) 

1 CONTINUE 
RETURN 
END 


FUNCT  ION  UZERO(U,X  ,'J PRIME, XPRIMF) 
r 

C*  44444444444444**.  ,444,44444* 

C r JNCTICN  U7ER0  COMPUTES  THr  VALUE  MU-ZERO  FOR  THE  TERM  * 

0 OMEGA  DOT  OMEGA  PRIMF  IN  CYLINDRICAL  COORDINATES.  * 

C*  * *444444444444444*  , 4444*4*44* 

C 

71  a 1.  - U**2 

72  = 1.  - U PR I ME 1 * 2 
OX  = X - XPRIME 

UZERC  = U *■  UPRIME  ♦ SQRT(Zl)  * S2RT(Z2)  * COS(DX) 

RETURN  S 

ENO  I 
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FUNCTION  SIG!1AG(9,Z,  AMUI,CHI  J,  AM JC,  CHID 


C 

E*  ******4**»***********»**»*** 

C Tu T S FUNCTION  DETERMINES  The  VALUE  0-  THE  “EVEN”  MACRO-  * 
C SCOPIC  SC  ATTrRI  NG  CROSS  S^TION.  * 

C*  **44**»**44*.»J»***|******»#4* 

n 

AMUNCT  = U7FR0(A‘lUItCHIJ,AMU»<,CHrL) 

PI  = 3.141G92bt>3E 
SIGMA  G = X3  (R,7,0>  /<»./pI 
L = LFCPPL ( IDUMMY) 

IF (L . LE.l ) RETURN 
CO  1 1=2, L, 2 

SIGMAG  = SIGMAG  ♦ ( 2* I +1 ) ♦ XS  (R , 7 , 1 ) * POLY ( I , AMUNO T > 
i n ./PI 
1 CONTINUE 
RETURN 
END 


FUNC7ICN  LFCPPL  ( T0J  ^Y) 

r, 

f**  *♦♦».»  4.  * 

C r JN  CT ICN  L EOF  DL  CALCULATES  THE  ORDF  . OF  THE  CROSS  SECTION* 
C rvr  ANSIOI  S . * 

<T*  * *4#»***4***4***»  t > f « I < 4 > * * » 

C 

COi'MC  M/EX  °A  f-'l)/  I°L 
LFOrr  L = I3L 
F ETUCN 
FNO 


FUNCT IOH  XS  <R»  Z» l > 

0 

f**  **»***.r*«,'.*44**  »>‘4-»»44,'»»* 

r -1‘JCTIOM  XS  CALCULATES  TH?  ANISOI  F;Or’lC  SCATTERING  CROSS  * 
C SECTION.  * 

/Ml  • 4444*4»»»»*«*»»  W»*»»  * »>  * 4 * 4 * 

c 

XS  = 1. 

TF(L.C-T.n)XS=0. 

FE7IJF  N 
FNC 


I 
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niMC'lCN  pCL  Y ( I » / ) 

C * a*..  .*,, 

C TWT9  Sl^:  CU7JV-  CALCUL  AT  FS  THE  I ’TH  L p > EN  DRF  POLYNOMIAL  * 

r fop  ru  ir  cut  \miue  o-  x.  * 

0“  «-•**  ♦■»♦*#**<,**♦*#,.»**»•,.,****, 

TOLY  = 1. 

J c ( I . EC  . C)  RETURN 
FCLY  = X 

IF(I.EO.l)  -.F  TURN 
CIM  = 1. 

PI  = X 
CO  1 L =2 » T 

p0i  Y = (.*?■  L-1)‘X‘  DT-(l.-l)*PIMl)/. 
cIf-’l  - pI 
PI  = POLY 
\ COI.T  I NUE 
TET'lf  |i 

END 


FUNC’  ION  SG  <R,"V  MM,  OHT) 

0 

r.  * ♦ *»*«*.><?».«  *.*»»«**♦* 

p r"  K!07  I OM  so  0 ETF RMIN rS  T«F  VALUE  OF  THE  ”FVEN"-PARTt  y * 
r E~>niCE  TFRM.  * 

C“  *•»•*♦>»»■»»*•*»•♦**♦«»♦****♦*♦# 

r 

SO  = 0. 

CFjUf  N 
FNO 


FUKOt  ION  S'J(R,  7,  6 MIJ,  P.H  T) 

C 

C*  ,.«*.»*.*»  «*♦*»*<►**.  > » » * i ««»«  + # 

o -UNCTION  EM  OCTCVtlNcS  THF  VALUE  OF  THE  "ODtV-P APITY  * 

n * 1 1 r r f tfpm.  * 

f 7*  »»^«*t»*»***»4**»»*»«*»*»*» 

c 

CU  r 0. 

FFlijr  N 

END 


FUNCTION  o*4DCOS(F|Zi  A MU » CHI ) 

0 

0*  ’ ***.*.t*»  >t»*»»4«*«»»*»  ****#• 

•“  r U NCT ICN  PNUCOS  CALCULATES  THF  BOUNTY  TFRHS  FOR  THE  FL-* 
p HfNT.  * 

p«  ♦ *¥•♦♦■***♦***»♦*  »*♦♦»**»** 

C 

rNOCOS  =0. 

RETURN 

FNO 
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FlIKC  ' ION  y I) 


c 

r*  •»«*+•»**  <•**•*»*■**•»*»«*«■♦♦♦* 

0 PIN'CTTC-N  P CALCUl  A rcr  THE  Alim  COORDINATE  CDr,k‘rSDOt;DING* 
C fTN-JIMFri?:  J‘|AL  XT  roor-'.OlNAT1'  Dr  THI-  "CANONICAL"  * 

C rL  -MFN1  . * 

V ♦ *♦■*■*  4 * » * • ’»•*♦♦•••«*»  4-  »»•*««■ 

r 

COFMf  N/M^o  Jpp.n/A  ,n,c,o 
F = A ♦ 3* XI 
FE1  <Jr  N 
FNO 


rUt-T-TICN  7(cTft) 

C, 

C*  * + * + *0***  + *»  + +.*  + xt**4l  + + ¥¥9 

C FUNCTION  z CfLCULA'EF  HE  ■»  COORDINATE  C0.W3P0U0ING  to  * 
o Tu  ~ NON-riHEUSIGML  FT  A ~OORQIMTF  n?  TH r •‘CANONICAL"  * 
C v '-mfnt  . » 

C*  **  + + **■*•*•  t*******  + #9*+  + + *.+ 

r 

CQKMF N/MAPCPAC/A, 7 ,c, d 
7 = r ♦ 0A£TA 
PETUF  N 
FNO 


rUNCT  T ON  AHJ  (XT) 

0-  * * » V ■»  -*  4 »*U*M**4»Ht  * + * t f.  + y * 9 

r 

c FUNCTION  A MU  CALCULATES  THF  mu  COOPOI NATF  CORRESPONDING  * 
C ft  THE  MFU-o:  MENSIONAL  XI  rOOPPTNA Tr  OF  THE  "CANONICAL"  * 
C Ft  EMFM  . * 

0 

<'0>'M f N/MA  PANGL  /F  »c  » 3,M 
AMU  = E ♦ F*XI 
RFlilF  N 
FNO 


FLU'C’  ION  CU(FTA) 

C 

****#*■•»»♦#»  **-»*♦» 

C FUNCTION  CNI  CALCULATES  THC  CNT  COORDINATE  COR1, FSPuNCING  * 
TO  "HE  NON-DIMENSIONAL  E’A  COOPOINATE  Dr  THE  "CANONICAL"  * 
0 cl  r MFNT  • * 

(*|  * »¥♦■«*«**•  >4  * * ***»J«**»* 


c 


OOHM<  M/MA  PANGL/r  ,c 
CH]  = G ♦ H " F T A 
PETUF  N 
FNO 


,G,H 
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* 


i 


SlT^f  UTII^  ‘U-’-’CLO,  AH-,AF>C,CLO,CHIf  CF^  C) 

* ♦■»  4 4 4 4 4 •»  f *4444  * - » 44k  4*4444 

TVJJS  Slir>r  OUTI'-r  COMPUT  PS  T“E  SCALING  PACTOPf  A FAC  AND  * 
~r"'  C 4b  ''TLL  AS  THE  “OfHrlIKATr  MAPPTg-,  \/ApIM_rS  US^O  TO  * 
■A‘ P f.COPC  IrJATF  J IN  Tur  l'|,CHI  ANCULAP  ° L A N F INTO  THP  XI,  * 
F’A  MCANf  HlCfcL"  PLAN r • THTS  POUTIN'  IS  USED  3v  FUNCTION  * 

r,njss2.  * 


ro-M(  n/ma  pang?/ a ,s  ,c , n 

A = (API  ♦ ALO'/?. 
n = ( A UI  - ALO)/?. 

C = (C^l  ♦ CLO)/?. 
r = (CHI  - ''LO)  / ?. 

APAC  = S 
CFAF  = 0 
RETlir  H 
END 


FUNCTION  CHI  2 ( ET  A ) 


4 4 4 4 4* 


44  *44444 


r r-.i motion  C.HI2  CALCJLATF;  TOE  chi  poo^inutf  CORRESPONDING* 
C TO  THI  FT  a COOPOINATF  O^  the  *TA  NO NT CAL'*  ELFMENT.  * 


4444444  44 


44*444*  4 


4 4 4 4 


cohhon/mapamg?/a»r»c» d 

CHI?  = C * n*  FTA 

CE7UF  N 

ENC 


FUNCTION  AMU? (XI) 


^UNCTION  A h u?  CALC  JLf ‘rcC  TwE  M'J  COO''!'! NA  T F CORRESPONDING  * 
t-1  THE  Nf  M-OlMrMGTONAL  XI  POOPPlNATP  0=-  T Hr  "CANONICAL”  • 
F,  FMFNT.  * 


C.OHMI  n/napang2/a,i 
AMU?  = A ♦ GVX I 
RETURN 
FNO 
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.■>  o o o .1  o o o r>  ->  o r>  n n o .■>  r>  „->  o o o o t o o •■>  t > o n «.">  r>  .*>  .*>  • ■>  o o 


FUNC1  I OK  3AMSS?(FCTl,FGT?,p,z,AMJ,CHI) 

( • ; 

»**<’*<»***»**»»,*t»»»*»«»»* 

F.'USS?  IS  A GA’JSS  -LEGF»ORF  MULTIPLE  r NTEG  RAT  I ON  ROUTINE  * 

•I*' T CH  USFS  AN  I POINT  Til  A OP.#'  TURl  RULr  0VEP  THE  MU,  CHI  * 

V*- ojARLC0  . • 

Jr (F, JtElAL)  IS  HF  GLOBAL  ANGLE  NOOr  MATRIX.  * f 

A r A C Ann  C F AC  ARE  COORDINATE  SCALING  rACT  ORS  SUPPLIED  BY  * 

SUBROUTINE  MAP?.  * } 

FUNCTION  GAULS?  USES  -JNCT70NS  AMU?, 0 HI?, 0,W, ^CTl ,FCT?»  ♦ ' 

Nr  , KE,  At  D LNN.  * 

COMMCN/HiLOW/ALDt  SO)  , A"I(FP)  ,Ci_n(£Q)  , CH1HK50) 

COM Ml N/SPACF/X IK, FT  AL/I  N OOF / LOCA . /O J A OPT  3 / 1 PO INT 

CO MM CN/NOTE S/M, N/rLEMENT/E/ANGFLM/NANG EL  ; 

COMMCN/TPI^LS/IT-riAL,  J” RIAL/ JE MAT  / JE(S0,16) 

INTEGER  E 

GAUSS?  = 0.  V 

n'TEPMINF  THf-  ‘.LORAL  ANGLE  NODE  OF  T-t  £ LOCAL  NOOF  JTRIAL 
U‘FD  IN  THE  CUTE*  ANGULAR  FLEMFNT  LOD». 

J = JF(F, JTRIAL) 

n“  TEPMINF  THE  NUMBER  OF  ANGULAR  ELEMrNT3  WHICH  CONTAIN 
Gl  09AL  NODE  J. 

NEL  = NE ( J , N,N ANGEL ,UF ) 

Ln  OP  OVFF  THF  ANGULAR  PLrurNxs  CONTAINING  GLOBAL  NODE  J 
AND  CALCULATE  THE  CO  ^ T =>  I BU  T I ON  EACH  M A <S  S TO  THE  INTEGRAL. 

DO  ? K=t, NEL 

DETFf  MINE  j HE  ELEMENT  NUMBER  OF  THE  K *TM  ANGULAR  ELEM. 

J *r  ^ 

KEL  = KE( J,K,N, NANGSL , IE> 

DETEf  MINE  THE  LOCAL  NONE  NUMBER  Dr  THE  GLOBAL  NODE  J 
IN  THE  K'TH  ANSU.AR  FLEMENT  CONTAINING  IT. 

>' 

LOCAL  = LNN  (J,KEL,N,JF) 

THE  VALUE  LOCAL  MUST  RF  PASSF.O  TO  FUNCTION  TP.IAL  ? ( SEE 
COMMON/LNOOE) . 

SET  L'P  Tmg  CANONICAL  MAPPING  OF  THE  < *TH  ANGULAR 
FLEMFNT  TO  BF.  INI  EGRATFO. 


5^ 


. I 


o c'  c c c t'  c ' t'  c;  c:  c:  o u o C' 


CAL  L MAP2<  ALO(  KE.)  , A-II  (KFL)  , AFAG,  :i_0  ( KFL  ) , GH  JHI  ( KEL)  , 
3 CF  AC ) 

FAC  = AFAC  • CFAC 


CAf-RY  OUT  NUMERICAL  OLK’DRITURE  OF  THF  K * TH  ANGULAR  EL. 

FAUST  = 0. 

00  1 J =1 » I cOIN  T 
XI  = 0(1) 

AHUI  = AN'JF(XI) 

00  1 J=  1 » I POINT 
ETA  = 0(J) 
rHIJ  - CHI^(FTA) 

F> A U S f = GAUSS  ♦ H(I)'  wrj)*ecT1(R,xfAMU,rHl>Af1IJIfCHIJ) 

5 *rCT2(XK,FTA.,XItFTA)*FAC 
1 CONTINUE 

GAUSS  ? = GAUSS  ? + GAUSS 
?.  CONTINUE 

F.ETUEN 
r NO 

F'JNCT  ION  Nc  ( J,  N,  .MANGEL  ,IE> 

'i  * » » « . . y 

~ ’TS  SUr;-'CU7lNE  CALCULATES  THE  UUH^-n  OF  ELEMENTS  CONTAIN* 
T" G GLrT'L  NJO--  J IN  THr  ANGULAR  00  5AIN.  ♦ 

J>  -*  « a • + .***  + ***  4 ***<•  *«*«»,  * * y ^ y 


CI-ENRION  I F ( L 0 » 1 6 ) 

1NTFCES  *-• 

I = 0 

GO  ? F = 1 , N ANGEL 
DO  1 k=i,«j 

tf(j.fo.:f(e,'<))  t=t*i 

T r ( I . f n ) GO  TO  ? 

1 FCMJNUE 
? CONTI  Mir 
* ME  = I 
RE  fljf  N 
FNG 

i 

» 

FUNCTION  LNN < J , K EL , N,  I c ) - 

*•»*♦*  V4  V«»>'4»44-*..I'  > *»■**»  » * * * * 

TuTS  SUOE  OUT  INF  CALCULATFS  THF  LOCAL  NO^F  NUMSFR  OF  GLOO-* 

Al  NODE  J IN  EL  EM  r NT  * 


niHEh'SION  I E ( C Of  16 ) 

00  1 1=1, N 

IF(IF (KEL, I) .EO. J) GO  TO  2 
1 CONTINUE 
? INN  = I 
RETUFN 

END  55 


k 

f. 

I 


• 4-  * i f ^ 


rum  3 (N  ( J ><  * N'XUGFL  >IE  > 

#44*  * * ♦ * 4 

r>  itlTc  SLr<*  CU’  1 *!-  f*\L»:UL#TES  T HF  WJM1P*  ^ 

^ l.  E L r r N ' rflMlMM>'3  L 0 ' A L f.'uor  I. 

• y 

nirr*  si?;  i f (‘>o*i: ) 

1 K'l  FT  TF  r 
K*  = C 

"0  ? r =1  ,i  vjr^'L 
f<3  1 1-1,'i 

IF  ( J.  c'1. ; ( r , I ) ) *' K = '» K ♦ H 

TF  ( . rn . <)  GL  TO  T 
1 c;OKT]Ni|r- 
->  r.Oi.1TU«JF 
f 1/  f r.  ? 

r.citjr  n 


=•  THf  <*TH  AU3U 

* •*  * t <■  * * 


Appendix  B 


Derivation  o 1’  I. lie  Weak  Korin  of  the  Transport  Equal. ion 


From  Chapter  II  we  have  the  functional 

CD 


r r w ^ 

F(u)  - JD[<n!l  V“',Ku(^Va)>  4 <u,6j(t> 


-2<*  ^M<«S„>-2<®S)>]^  +jfj/|£-aM|  «* 
» <*,S>  * f -Hi)  $£)*■& 


JiJL 


<ir 


(23) 


where 


and  * stands  for  the  complex  conjugate. 


We  want  to  show  that  F(u)>FM  where  U = V<4*J  andvj^O.  From  Eq  (23) 
a^>ve,  we  have 


(B-1 ) 


(B-2) 


+ f <(^'V^),Ku(iL'V>i)><i£ 

D 

+ f <(£  ^*l)  > Ku.(4'^’V;)>4  (B-3) 

+ /D<(^k ' vn)  > Ku(i  yl)>4 

*ID<W>i* 

+ /„<*.  Gj1>4* 

+ /n<'1,G,t>fr 


(B-4) 


(B-5) 


(B-6) 


(B-7) 
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+ J D<fl'G)7\> 

(B-8) 

-2  [ <. a ^,Kusa>  U 

(B-9) 

-2JD<Z<”l,KuS^>h 

(B-10) 

-2jD<*,s3yh 

(B-11 ) 

_2/d<’ >'s>>ir 

(B-12) 

+ JRfji  U-»«|  i-x. 

(B-13) 

(B-l4) 

+ / \-ti  | ^ <Lt 

JR 

(B-15) 

It  can  be  seen  that  Eq  (B-l)  + Eq  (B-5)  + Eq  (B-9)  + Eq  (B-ll)  + Eq 
(B-13)  equal  fM.  Since  the  operators  G^,  G^,  K^,  and  are  self- 
ad  joint  , that  is, 

<m,  Sj  3®>  = <g,*(4,  ,®> 


- <3W,G,4«>* 

< 

we  have  terms  j 

t 

(B-2)  = (B-3)* 

and  , 

(B-6)  , (B-7)*  * 

1 
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>— a . /»  ■ 


. ' Eg 


I 


Therefore,  terms 


(b-,/)  + (B-3)  = (U-3S ) * + (B-3)  = P.RcCb-3) 


and  we  find  that 


F(%n)  * f(^) 


+ 2 / [Re<-8,-^yl . Ktt(JL-V^)>  + Re<>!  ,G. 
LJD  (B-  3)  (B - 7) 


- < Jl  V>l,  fVS*>-  <*i,S,>]  JLr 
(B-lt>)  ( B - 12 ) 

+ ^d/J-'21WI  ^ 

R - (B'  14)  J 

+ /Dt  <'l,6,v|>]Ar 


(B-4) 


(B-8) 


+ I Yl't  ir 
K - (r-im 


" * (B~  15)  J 

Because  the  operators  G , G , K , and  K are  all  positive-definite,  the 

g u g u 

terms  (B-h),  (B-8),  and  (B-15)  in  the  bracketed  (C)  term  are  positive 
as  long  as  Vf  is  not  equal  to  zero.  Thus,  the  (C)  term  is  positive. 

Similarly,  the  operators  G^,  G^,  K^,  and  are  real  operators; 
they  give  real  results  when  operating  on  real  functions.  So,  as  long 
as  we  use  real  functions  *r  and  V|  , we  can  drop  the  two  Re's  in  the 
bracketed  term  (B). 

Using  the  property 
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I 


U J Q 


+ <£  <(i' *(>£})*,  ,-f>lr  CB-16) 
R 


the  (b-3)  and  (B-10)  terms  in  the  bracketed  term  (B)  can  be  written  as 


/ <n,--a-'?|'K^k-v  <(']]  > Ai 


+JR  < .KJa-v  Vj>  ir  (b-3)  • 


/<>),-£' V(K„  SJ>iv 


+ /R  <(**)  VI,  K„s*>^ 


(B-10) 


Bracketed  term  (B)  thus  becomes 


2/  <>1  -41  v[Ka(^.^  <//)]>+  <-n,G3  (^> 

r\ 


+ <>1  v(kuSJ>- <^,S3>  J 1* 


-<(•£**)  v],k*Su> 

. 


i 


I 


~2Jd  [<n^MK*(^'^)] + G3  ^ [Sj  ~ ^-■57(Kusu}  j>] 

+ 2^4  vl[|-a-s|  ^ * fe  !)[Ku(a-p'/')-KuSu]j  Jaiv 

Eqs  (l/),  (l8a),  ajid  (l8b)  of  Chapter  II  are  now  applied  to  force  term 
(B)  to  zero.  Thus,  F(4/+y|)  is  indeed  greater  than  F(^) . The  resulting 
weak  form  of  the  transport  equation  corresponding  to  its  variational 
formulation,  Eq  (2J),  is 

JD  l<k'Vv|,Xj4.-v^)>+  <>1iGj4'>]  <U 

= [<C-S=-^7),KUSU>  + <‘*/l,S3>]  cir  (2*0 
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